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ABSTRACT 


An  HF  direction-finding  (DF)  system  employing  sampled-aperture  techniques  on  high- 
latitude  signals  is  simulated.  The  propagation  is  treated  as  scatter  from  moving  patches  of  field- 
aligned  F-region  irregularities  modelled  as  vertical  “thin-wire”  scatterers  in  a  horizontally 
stratified  ionosphere,  plus  a  weak  sporadic-E  (Eg)  reflection  from  a  horizontal  mirror-like  surface. 
The  resultant  simulated  received  signals  are  analysed  by  plane  wave  fitting  and  multiple-signal 
classification  (MUSIC)  DF  algorithms. 

Simulations  of  DF  operations  on  scattered  signals,  using  system  parameters  similar  to  those 
used  for  actual  high-latitude  observations,  reveal  results  similar  to  those  observed.  The  observed 
directions  are  an  indication  of  the  locations  of  the  scattering  patches  rather  than  the  transmitter. 
The  detection  of  a  patch  depends  strongly  on  its  location,  the  radio  frequency,  and  the  ionospheric 
parameters.  When  a  patch  is  seen,  the  time  history  of  estimated  directions  indicates  its  width,  as 
well  as  its  location.  MUSIC  provides  a  picture  of  the  width  and  number  of  scattering  patches 
present  in  a  short  time.  Plane-wave  fitting  tends  to  give  only  the  directions  for  the  patch  having 
the  strongest  scattered  signal  and  requires  a  longer  time  for  the  successive  direction  estimates 
(spread  over  the  angle  subtended  by  the  patch)  to  indicate  the  patch  width.  When  too  many  signal 
directions  are  assumed,  MUSIC  gives  false  directions  which  correspond  to  the  steered  directions 
of  the  array  for  which  the  larger  sidelobes  are  pointed  toward  the  strongest  signals.  When  a 
weaker  Eg  signal  is  present  with  the  scattered  signal,  MUSIC  can  see  it  clearly  and  provide  a  good 
direction  estimate,  provided  array  dimensions  are  not  so  large  that  narrow  beamwidths  coupled 
with  large-sidelobe  effects  prevent  the  weaker  Eg  signal  from  being  seen. 
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RESUME 


Un  dispositif  radiogoniom^trique  haute-fr6quence  k  6chantillonage  d’ouverture  est  simule 
pour  le  cas  des  signaux  que  Ton  retrouvre  aux  haute  latitudes.  Dans  la  simulation  on  considere 
que  la  propagation  des  signaux  est  le  resultat  del  la  dispersion  des  signaux  par  des  zones  non- 
uniformes  et  mouvantes  de  multiples  conducteurs  minces  dans  la  couche  F  alignes  sur  le  champ 
magn^tique  terrestre.  Ces  zones  sont  representes  dans  la  simulation  par  des  fils  conducteurs 
minces  agissant  comme  elements  dispersifs  dans  une  ionosphere  stratifiee  horizontalement, 
auxquels  on  ajoute  I’effet  d’une  reflection  faible  provenant  de  la  couche  sporadique  E  (Eg) 
represent^  dans  la  simulation  par  un  miroir  horizontal.  Les  signaux  resultants  sont  analyses  par  le 
moyen  d’algorithmes  impliquant  ajustement  d’ondes  planes  et  classification  de  signaux  multiples 
(appeles  MUSIC  pour  Multiple  Signal  Classification). 

Les  resultats  obtenus  a  la  suite  de  simulations  de  systemes  radiogoniometriques  operant  sur 
signaux  de  dispersion  permettent  d’obtenir  des  resultats  conformes  aux  mesures  aux  hautes 
latitudes.  Les  resultats  del  la  goniometrie  revelent  que  les  directions  de  provenance  des  signaux 
correspondent  aux  zones  de  dispersion  at  non  a  la  position  de  I’emetteur.  La  detection  d  une  zone 
de  dispersion  depend  principalement  de  sa  position,  de  la  frequence  radio  utilisee  et  des 
parametres  physiques  de  I’ionosphere.  L’historique  de  I’observation  de  la  direction  d  une  zone  de 
dispersion  permet  d’en  deduire  la  largeur  et  la  position.  Le  systeme  MUSIC  permet  d  obtenir  une 
representation,  pour  une  courte  p^riode  de  temps,  de  la  largeur  et  du  nombre  de  zones  de 
dispersion.  La  technique  d’ajustment  d’onde  plane  permet  seulement  de  determiner  les  directions 
successives  de  la  zone  dispersive  produisant  la  plus  grande  intensity  de  signal.  En  outre,  cette 
technique  exige  un  temps  plus  considerable  de  calcul  en  ce  qui  a  trait  k  I’estimation  (s  etalant  sur 
Tangle  sous-tendu  par  la  zone  de  dispersion)  des  directions  successives  des  signaux  provenant  des 
zones  de  dispersion  pour  en  arriver  au  calcul  de  la  largeur  de  la  zone.  Quand  trop  de  directions  de 
signaux  incidents  sont  presuppos^es  le  systeme  MUSIC  produit  de  faux  resultats.  Les  directions 
ainsi  obtenues  correspondent  Talignement  de  lobes  lateraux  avec  les  signaux  les  plus  intenses. 
Le  systeme  MUSIC  permet  la  determination  de  la  direction  de  signaux  faibles  Eg  (parmi  les 
signaux  forts  de  dispersion)  pouvu  que  le  dimensions  du  reseau  d  antennes  soient  telles  que  le 
faible  etalement  du  faisceau  principal  combine  aux  effets  dus  aux  lobes  lateraux  n’en  previennent 
la  detection. 


EXECUTIVE  SUMMARY 


The  operation  of  a  sampled-aperture  HF  direction-finding  system  at  high  latitudes  is 
considered.  The  irregular,  drifting  ionosphere  present  at  high  latitudes  provides  the  conditions  for 
the  creation  of  patches  of  field-aligned  irregularities,  which  scatter  radio  waves.  This  scattering 
can  cause  radio  signals  to  be  received  from  directions  well  removed  from  the  great-circle  bearing 
of  the  transmitter.  A  feature  which  may  also  be  present  some  of  the  time  that  scattering  occurs  is 
sporadic-E.  Sporadic-E  reflections  result  in  radio  propagation  close  to  the  great-circle,  allowing 
an  accurate  bearing  estimate  to  be  made. 

A  model  is  developed  for  the  scatter-supported  ionospheric  radio  propagation,  along  with 
the  sporadic-E  propagation.  The  major  assumptions  underlying  the  model  are  the  concept  of 
‘scattering  centers’,  or  single  locations  where  all  the  scattering  for  a  particular  path  is  assumed  to 
take  place,  and  a  horizontally  stratified  ionosphere.  Although  these  assumptions  are  an 
oversimplification  of  the  actual  high-latitude  situation,  they  do  permit  a  first-order  analytic  model 
to  be  developed,  from  which  some  insight  into  the  actual  situation  can  be  gained.  The  scattering 
centers  are  assumed  to  exist  in  ‘patches’  of  limited  horizontal  extent,  consistent  with  knowledge 
of  the  high-latitude  ionosphere. 

The  direction-finding  operation  is  also  modelled.  This  includes  a  two-dimensional  array  of 
antennas,  and  a  direction-finding  algorithm  which  obtains  direction  estimates  from  the  signal 
samples  received  at  the  array  elements.  Two  algorithms  are  modelled:  a  plane-wave  fitting 
algorithm,  which  obtains  a  single  direction  corresponding  to  the  best  plane-wave  fit  to  the 
received  signal  samples,  and  a  MUSIC  algorithm,  which  is  capable  of  obtaining  a  number  of 
signal  directions  from  a  set  of  signal  samples. 

DF  system  operation  is  simulated  by  running  the  model  a  number  of  times,  with  the 
scattering  patch  locations  shifted  a  small  amount  each  time,  to  simulate  ionospheric  drift  motions. 
The  resultant  time  histories  of  estimated  bearing  and  elevation  are  then  examined  to  determine  the 
appearance  of  the  scattering  and  sporadic-E  features,  and  the  performance  of  the  system. 

A  particular  set  of  simulations  were  used  to  study  the  time  history  of  the  observed  signal 
directions  from  scattering,  for  various  scattering-patch  geometries  and  propagation  parameters 
such  as  ionospheric  electron  density.  The  transmitter  and  receiver  conditions  simulated  were 
Rimilar  to  those  of  the  Kestrel  measurements  conducted  at  Alert  for  transmissions  from  Thule,  in 
November  1990.  The  Kestrel  measurements  used  a  plane-wave  fitting  technique  to  obtain 
direction  time  histories.  Likewise,  the  simulations  used  plane-wave  fitting  to  derive  direction  time 
histories  which  may  be  compared  to  those  of  Kestrel. 

It  was  found  from  these  simulations  that  the  signal  from  an  extended  patch  fits  a  plane  wave 
quite  well  most  of  the  time,  with  the  arrival  direction  coming  from  anywhere  inside  the  patch. 
Consecutive  estimates  yield  different  directions  that  are  distributed  over  the  angle  subtended  by 
the  patch  at  the  receiving  array,  so  that  a  rough  estimate  can  be  obtained  of  the  size  of  the 
scattering  patch.  The  appearance  of  the  time  histories  was  similar  to  many  of  the  Kestrel  records. 
The  spread  in  bearing  observed  in  the  simulations  is  substantially  greater  than  the  spread  in 
elevation,  a  feature  noted  in  many  of  the  Kestrel  records  as  well.  The  ionospheric  propagation 
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determines  the  portion  of  the  ionosphere  wherein  a  scattering-center  can  be  located  which  will 
contribute  to  the  propagation.  For  a  horizontally  stratified  ionosphere,  the  contributing  scattering 
centers  are  limited  to  a  strip,  perpendicular  to  the  direct  transmit-receive  path  and  center^  on  its 
midpoint.  For  an  irregular  ionosphere,  the  allowed  contributing  region  could  be  quite  different. 
With  moving  scattering  patches,  the  patches  are  seen  only  where  they  pass  through  the  allowed 
region.  When  two  or  more  patches  are  present,  the  patch  which  has  the  strongest  composite 
scattered  signal  at  the  receive  array  will  control  the  estimated  direction.  W^th  random  scattering 
processes,  the  signal  received  from  any  one  patch  can  vary  substantially  from  estimate  to 
estimate,  so  that  when  two  patches  are  fairly  close  in  their  mean  signal  power,  the  estimated 
direction  jumps  back  and  forth  between  patches,  allowing  both  of  them  to  be  seen. 

When  the  plane-wave  fitting  procedure  was  replaced  by  the  MUSIC  algorithm,  the 
scattering  simulations  showed  a  generally  improved  result.  Because  a  number  of  signal  directions 
are  obtained  for  a  single  estimate,  it  becomes  possible  to  get  a  clear  picture  of  the  number  and 
extent  of  the  scattering  patches  with  only  a  very  short  sequence  of  estimates.  One  problem  that 
was  observed  with  MUSIC,  however,  is  the  occasional  presence  of  false  directions  that  occur 
when  too  many  assumed  signal  directions  are  used.  These  false  directions  correspond  to  array 
directions  for  which  a  major  sidelobe  of  the  array  is  directed  toward  the  stronger  signals. 

The  simulations  which  include  sporadic-E  as  well  as  scattering  propagation  clearly  show 
that  the  MUSIC  algorithm  can  be  used  quite  effectively  to  enable  relatively  weak  sporadic-E 
signals  to  be  seen  in  the  presence  of  strong  scattered  signals.  This  holds  some  promise  for  high- 
latitude  direction  finding,  since  the  sporadic-E  signals  are  expected  to  be  quite  representative  of 
the  true  bearing  to  the  transmitter  while  the  scattered  signals  are  not.  The  masking  of  weaker 
sporadic-E  signals  by  the  array  sidelobe  response  to  stronger  scattered  signals  is  a  possibility, 
especially  when  the  array  beamwidths  are  much  narrower  than  the  range  of  directions  the 
scattered  signals  arrive  from.  It  is  therefore  important  to  avoid  large-aperture  sparsely-filled 
arrays  which  have  strong  sidelobes.  Three  array  geometries  were  tested:  a  the  Kestrel  array 
consisting  of  7  elements  with  a  constant  radius  of  25  m,  the  10-element  middle  portion  of  the 
(Vortex)  log-spiral  array  currently  installed  at  CRC  with  a  spiral  radius  varying  from  27.5  to  175.0 
m,  and  the  8-element  inner  portion  of  the  Vortex  array  with  a  spiral  radius  varying  from  14.8  to 
62.9  m.  The  best  simulated  result  were  attained  for  the  inner  8  elements  of  the  Vortex  array.  With 
this  array,  at  an  operating  frequency  of  8  MHz,  it  was  possible  to  obtain  accurate  directions  for 
sporadic-E  signals  21  dB  below  the  scattered  signals.  A  number  of  the  assumed  signal  directions 
are  taken  up  by  the  signals  scattered  by  the  extended  scattering  patches,  so  it  is  important  to  have 
enough  assumed  signal  directions.  Too  many  directions  will  result  in  false  signal  direction 
estimates,  as  a  result  of  the  array  sidelobe  effect.  Since  the  number  and  extent  of  the  scattering 
patches  are  not  known  beforehand,  it  is  important  that  an  operational  system  try  various  numbers 
of  assumed  signal  directions,  in  order  to  obtain  a  good  result. 
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1.0  INTRODUCTION 


The  simulation  studies  presented  here  are  an  attempt  to  model  the  operation  of  a  sampled- 
aperture  HF  direction-finding  (DF)  system  operating  at  high  latitudes.  The  goal  of  this  work  is 
two-fold:  to  obtain  propagation  models  that  are  representative  of  high-latitude  conditions,  and  to 
test  various  DF  methods  against  those  models  to  establish  their  utility. 

HF  direction-finding  systems  depend  on  radio  propagation  and  therefore  the  condition  of  the 
ionosphere  through  which  HF  signals  must  pass  in  order  to  be  seen.  At  high  latitudes,  the 
ionosphere  is  highly  changeable,  so  that  its  effect  on  HF  radio  propagation  is  especially 
pronounced. 

The  high-latitude  ionosphere,  which  includes  both  the  polar  cap  and  auroral  regions,  is 
governed  by  the  interaction  between  the  solar  wind  and  the  earth’s  magnetosphere.  The 
ionosphere  is  located  above  the  neutral  atmosphere.  The  earth’s  magnetic  field  lines  pass  through 
the  ionosphere,  and  extend  upward  into  the  magnetosphere.  Any  interaction-generated  motions 
which  occur  in  those  field  lines  are  mapped  down  onto  the  ionosphere.  These  motions  are  thought 
to  take  the  form  of  a  two-cell  convective  pattern,  whose  direction  is  mainly  antisunward  across 
the  polar  cap  region  covering  the  poles,  and  sunward  at  the  auroral  and  subauroral  regions 
equatorward  of  the  polar  cap  [1,  2].  The  ionosphere  in  the  polar  cap  tends  to  be  more  irregular 
during  darkness,  since  the  most  steady  source  (photoelectron  production)  and  the  main 
irregularity  decay  mechanism  (E-layer  conductivity)  have  ceased  at  those  times.  A  prominent 
feature  of  the  nightside  polar  cap  ionosphere  is  the  large-scale  blobs  of  enhanced  electron  density 
convecting  antisunward  from  the  dayside  [3].  At  the  same  time,  the  moving  gradients  associated 
with  these  blobs  give  rise  to  short-lived  small-scale  irregularities  whose  locations  move  with  the 
same  convective  flow  [3]. 
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Figure  1.  Electron-density  contour  plot  of  the  high-latitude  ionosphere,  obtained  by  the 
Chatanika  incoherent  radar,  for  Feb.  20, 1981,  during  local  night  (1130-1144  UT)  [4]. 
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Figure  1  is  a  contour  plot  of  the  high-latitude  ionospheric  electron  density,  taken  by  the 
Chatanika,  Alaska  incoherent  radar  facility  [4].  The  figure  is  a  vertical  north-to-south  cut  through 
the  night  portion  of  the  ionosphere,  showing  the  subauroral  region  (61  to  63.5°  N),  the  auroral 
region  (63.5  to  66.5°  N),  and  a  portion  of  the  polar  cap  (above  66.5°  N).  The  strongest  horizontal 
gradients  are  seen  in  the  auroral  region.  Significant  large-scale  irregular  features  are  also  seen  in 
die  polar  cap  region  which  is  of  major  interest  to  the  present  work. 

The  implications  for  HF  radio  propagation  are  substantial.  Radio  waves  will  be  reflected 
obliquely  from  large  and  medium-scale  irregularities  of  several  km  extent  or  greater  (one  Fresnel 
zone  for  HF  radio  waves).  They  also  will  be  scattered  from  smaller-scale  irregularities  which  have 
sizes  less  than  one  km.  The  larger-scale  features  are  expected  to  be  sufficiently  long-lived  that  in 
their  lifetimes  they  will  move  substantially  with  the  general  high-latitude  convective  flow[5].  The 
smaller  irregularities  last  several  seconds  or  less[6].  These  irregularities  exist  in  patches  in  the 
vicinity  of  the  large-scale  gradients  which  drift  with  the  convective  flow. 

Two  complementary  models  for  high-latitude  radio  propagation  have  been  developed  in  the 
current  study  program.  The  first  model  is  based  upon  reflection  of  radio  waves  from  large  and 
medium  scale  irregularities  and  treats  the  irregular  ionosphere  as  a  bumpy-mirror  reflecting 
surface.  The  bumps  were  modelled  as  fixed  and  occurring  in  patches  which  moved  with  a  constant 
horizontal  motion  typical  of  high-latitude  drift  motions.  The  second  model  is  based  upon  aspect  - 
sensitive  scattering  from  the  small-scale  irregularities.  These  models  have  been  presented  in 
technical  memoranda  [7,  8]. 

A  currently  available  source  of  high-latitude  DF  data  is  that  of  the  Kestrel  program  [9, 10]. 
This  program  consisted  of  a  two-week  measurement  campaign  utilizing  seven  elements  of  the 
Alert  pusher  array  together  with  an  Andrew  Corp.  sampled-aperture  receiving  and  processing 
system.  This  system  performed  a  real-time  phase-front  analysis  on  the  sampled  signals  from  the 
array,  to  obtain  the  best-fit  plane-wave  bearing  and  elevation  which  were  recorded  along  with  a 
measurement  quality  index  denoted  as  the  WSPI,  several  times  a  second.  These  recordings  form  a 
permanent  record  against  which  the  predictions  of  models  can  be  tested.  The  most  helpful 
measurements  for  establishing  high-latitude  propagation  conditions  are  those  taken  on  the  8.05 
MHz  transmissions  from  Thule,  667  km  distant.  They  provide  a  propagation  path  lying  well 
within  the  polar  cap  at  all  times  of  day. 

Both  the  bumpy-mirror  and  aspect-sensitive  scattering  models  were  tested  against  a  plane- 
wave  fitting  DF  procedure,  using  an  array  and  transmitter  geometry  similar  to  that  of  the  Kestrel 
data  for  Thule  [7,  8].  This  provided  a  means  of  comparing  the  predictions  of  these  models  against 
observed  data.  Both  models  predict  a  spread  in  the  consecutive  arrival  directions  found  with  a 
plane  wave  estimator,  and  a  change  in  direction  corresponding  to  the  general  motion  of  a  bumpy 
region  or  region  of  scatterers.  These  predictions  are  consistent  with  the  Kestrel  observations. 
However,  the  bumpy-mirror  model  by  its  nature  imposes  a  large  spread  in  elevation  when  the 
spread  in  azimuth  is  large  [7],  while  the  scattering  model  imposes  a  much  smaller  elevation 
spread  because  of  the  selective  nature  of  aspect  sensitivity  [8].  In  this  regard,  the  Kestrel  data  is 
more  supportive  of  the  scattering  model  [9].  Also,  the  scattering  model  predicts  that  different 
scattering  regions  will  be  visible  in  a  sequence  of  consecutive  direction  estimates  [8],  which 
appears  to  be  the  case  at  least  part  of  the  time.  The  bumpy-mirror  model  predicts  a  more  abrupt 
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transition  in  the  visibility  of  two  moving  disconnected  bumpy  regions  as  their  relative  strengths 
change  [7].  Nonetheless  the  bumpy-mirror  model  does  predict  one  unique  feature  which  is  seen  in 
the  data  at  times:  the  occasion^  clustering  of  directions  (associated  with  randomly-occurring 
strong  focusing  from  a  single  reflecting  location).  The  bumpy-mirror  model  is  not  considered 
further  in  this  report. 

One  propagation  mode  not  included  in  the  models  [7,  8]  but  seen  in  the  Kestrel  data  is  that 
of  sporadic-E  (Eg)  [9].  This  mode  occurred  only  rarely  in  this  data,  and  lasted  only  for  several 
minutes  at  a  time.  However,  the  direction  estimates  obtained  with  this  mode  were  much  more 
accurate  than  the  estimates  obtained  at  other  times.  The  DF  plane-wave  fitting  technique 
employed  in  this  data  permitted  the  Eg  mode  to  be  seen  only  if  its  signal  strength  was  comparable 
to  or  exceeded  that  of  other  modes.  In  view  of  the  improved  accuracy  for  this  mode,  other 
sampled-aperture  DF  techniques  should  be  considered  if  they  allow  the  Eg  mode  to  be  seen  when 
it  is  lower  in  strength  than  the  remaining  modes.  A  class  of  algorithms  typified  by  MUSIC  [11] 
fall  into  this  class.  These  algorithms  permit  a  number  of  arrival  directions  to  be  ascertained  from  a 
single  set  of  data,  even  when  the  corresponding  signal  levels  are  quite  different  in  their  strengths. 


17-H0U-90 
Start  Titia  1900 


AZIMUTH  &  ELEVATION  SEQUENCES 


LIMITS <az  0-3S9  uspi  10> 


80S(k.^  kHz 


Figure  2.  Kestral  azimuth  and  elevation  time  histories  for  the  Thule  8.05  MHz  transmissions. 


An  example  of  the  Kestral  data  is  given  in  Figure  2,  which  shows  the  tracking  of  a  moving 
scattering  region,  the  spread  of  its  estimated  arrival  directions,  and  the  brief  appearance  of  the  Eg 
mode  from  42  to  50  minutes  from  the  start.  The  Eg  mode  is  characterized  by  its  closeness  in 
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azimuth  to  the  great-circle  direction  (194°  E  of  N)  and  its  lower  elevation  angles  (as  compared  to 
the  F-layer  returns  seen  at  other  times). 

The  present  formal  report  includes  a  description  of  the  scattering  model  used.  The  model  is 
extended  beyond  the  single-sample  model  presented  previously  [8]  so  that  it  may  be  used  with  DF 
algorithms  such  as  MUSIC  which  require  more  than  one  sample  of  the  signals  across  the  antenna 
array. 


The  predictions  of  the  extended  scattering  model  for  sampled-aperture  DF  arrays  are 
presented,  using  both  a  plane-wave  fitting  DF  procedure  and  the  MUSIC  algorithm.  These 
predictions  are  of  prime  interest,  in  determining  the  validity  of  the  scattering  model  in 
comparisons  with  acquired  data  and  in  assessing  its  effect  on  the  performance  of  the  two  DF 
algorithms. 

In  a  second  series  of  simulations  also  included  in  this  report,  a  sporadic-E  propagation  mode 
has  been  added  to  the  scattering.  The  prime  concern  here  is  with  the  performance  of  a  DF  system 
in  the  presence  of  an  ‘extended-source’  or  ‘many-path’  propagation  mode  such  as  scattering, 
along  with  the  fixed  single-path  Eg  mode.  The  MUSIC  algorithm  is  tested  for  its  ability  to  see  the 
Eg  mode  when  its  signal  is  much  weaker  than  that  of  the  scattered  signal.  Various  algorithm  and 
antenna  configurations  are  tested  to  determine  which  are  most  effective  under  these  conditions. 
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2.0  THE  PROPAGATION  MODEL 


2.1  Overview 

The  propagation  model  is  comprised  of  two  parts:  F-layer  scattering  and  sporadic-E 
propagation.  A  flat  earth  is  assumed.  The  F-region  is  assumed  to  be  horizontally  stratified  and 
parabolic  in  its  height  profile.  The  effect  of  collisions  and  the  magnetic  field  on  radio  wave 
propagation  are  neglected,  resulting  in  a  real  isotropic  refractive  index  [12].  The  E-region  is 
assumed  to  be  a  flat  horizontal  reflector,  with  a  specified  power  reflection  coefficient.  A 
transmitting  and  a  receiving  location  are  specified,  along  with  a  receiving  array  geometry  and  the 
elevation  patterns  of  the  transmitting  and  receiving  antennas. 

The  direction-finding  system  modelled  is  a  sampled-aperture  array,  with  either  a  plane-wave 
fitting  algorithm  or  the  MUSIC  [11]  algorithm  applied  to  the  signal  samples.  The  plane-wave 
fitting  algorithm  requires  only  a  single  sample  from  each  antenna,  and  provides  a  single  direction 
estimate.  The  MUSIC  algorithm  requires  a  number  of  consecutive  samples  be  taken  from  each 
antenna,  and  permits  more  than  one  signal  direction  to  be  estimated  fi'om  these  samples. 

The  modelled  scattering  occurs  in  the  F-layer,  in  specified  regions  or  patches  which  contain 
a  specified  density  of  vertically-aligned  ‘scattering  centers’.  These  scattering  centers  are 
introduced  in  order  to  replace  the  multiple-scattering  problem  by  the  more  tractable  single¬ 
scattering  problem.  This  approach  is  partially  justified  by  noting  that  for  vertically-aligned 
elongated  scatterers  of  small  horizontal  dimension,  a  singly-scattered  wave  has  the  same  elevation 
angle  as  the  original  wave  and  a  completely  arbitrary  azimuth,  i.e.,  the  incident  wave  is  scattered 
equally  into  all  azimuths.  As  a  result,  multiple  scattering  also  preserves  the  elevation  angle  and 
the  energy  is  scattered  equally  into  all  azimuths,  regardless  of  the  number  of  scatters.  The 
scattering  centers  replace  the  scattering  that  occurs  at  several  places  along  the  path  of  a  particular 
ray  by  a  single  scattering  occurring  at  the  location  of  the  scattering  center.  The  propagation 
process  is  described  further  in  Section  2.3.  The  energy  intercepted  by  a  scattering  center  is 
scattered  equally  into  all  azimuths,  while  the  elevation  angle  is  preserved.  A  height-dependent 
‘scattering  width’  representing  the  scattering  efficiency  is  introduced  for  each  scattering  center. 
Given  a  particular  scattering  center,  under  the  assumption  of  a  horizontally  stratified  ionosphere 
the  potential  propagation  paths  have  a  common  projection  on  the  horizontal  plane,  consisting  of  a 
straight  line  from  the  transmitter  to  the  base  of  the  scattering  center  and  another  from  the 
scattering  center  to  the  receiver. 

Each  scattering  center  generated  by  the  model  is  considered  separately  in  terms  of  its 
contribution  to  the  total  signal  seen  at  the  receiving  array.  The  MUSIC  algorithm  requires  a  series 
of  consecutive  signal  samples  in  order  to  obtain  a  non-singular  covariance  matrix  from  which  the 
directions  are  estimated.  In  practice  these  samples  would  be  close  in  time,  perhaps  a  fraction  of  a 
millisecond  apart.  The  scatterers  themselves  may  have  time  constants  of  the  order  of  a  tenth  of  a 
second  or  more.  In  order  to  simulate  the  time  changes  and  generate  consecutive  samples,  the 
contribution  of  each  scattering  center  is  modulated  in  its  phase  and  amplitude  by  narrowband 
white  noise,  the  bandwidth  of  which  is  set  according  to  the  Doppler  spreading  anticipated  on 
high-latitude  HF  circuits.  The  actual  locations  of  scattering  centers  are  considered  constant  over  a 
series  of  consecutive  signal  samples. 
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The  simulation  is  performed  for  a  number  of  estimates,  to  simulate  direction-finding  on  a 
particular  signal  over  a  period  of  time.  These  estimates  are  taken  from  separate  sets  of  samples, 
between  which  the  patch  locations  have  shifted  by  a  specified  constant  amount  in  order  to 
simulate  high-latitude  ionospheric  convection.  For  each  set  of  samples,  the  actual  locations  of  the 
scattering  centers  within  each  patch  are  generated  anew.  The  parameters,  such  as  scattering-center 
density  and  patch  dimensions,  are  assumed  to  remain  constant  over  a  simulation  run. 


2.2  Generating  Scatterers 

In  the  flat-earth  geometry  of  the  model,  a  scattering  center  location  is  expressed  by  its 
horizontal  coordinates  x,y.  Each  center  is  vertically  oriented,  and  can  be  considered  to  exist  at  all 
heights  h  .  The  scattering  regions  or  patches  are  assumed  rectangular  so  that  they  can  be  specified 
by  their  x  and  y  coordinate  ranges  jcq  to  » Jo  to  • 

For  a  given  mean  density  of  scattering  centers  v,  the  probability  of  a  scatterer  existing  in  an 
infinitesimally  narrow  strip  dx  wide  between  yQ  and  is  given  by 


pix)dx  =  viyQ-yi)dx  .  (t) 

This  results  in  the  scattering  centers  being  distributed  in  their  x-coordinate  according  to  a  Poisson 
distribution,  so  that  the  distance  d^^  along  the  x-axis  from  jtq  to  the  first  scattering  center,  and 
between  succeeding  centers,  is  governed  by  the  distribution  function 

Pid>  d^)  =  exp  (-V  (yo - y^)  d^)  (2) 


The  Monte  Carlo  method  is  used  to  compute  the  distance  along  the  x-axis  to  the  first  scattering 
center  and  between  succeeding  centers.  Equation  (2)  is  inverted,  and  P  {dt.  d^)  is  replaced  by 
independent  computer- generated  random  numbers  uniformly  distributed  between  0  and 

1.  The  x-axis  separations  are  then  given  by 


d^i) 


vCyo-Ji) 


(3) 


where  the  index  /  indicates  the  ith  scattering  center.  Thus  the  x-coordinate  of  the  ith  center  is 
given  by 

^,(0  =  +  for/=l 

(4) 

x^(i-l)  +d^ii)  for/>l 

The  x-coordinates  are  generated  sequentially  according  to  equation  (4),  until  x^  is  exceeded.  The 
corresponding  y-coordinates  are  distributed  uniformly  between  yQ  and  y^  so  that 

ys(0  =  yo+  (yo-yo^rand^^) 

where  the  N'^and^^^  ^  likewise  independent  computer-generated  random  numbers  between  0 
and  1. 
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2.3  Propagation  Geometry 


The  propagation  geometry  of  the  model  is  illustrated  in  Figure  3.  For  a  particular  successful 


Figure  3.  Propagation 
geometry  showing  transmitter, 
scattering  center,  and  receiver 
locations. 


vertical 


path,  the  ray  leaves  the  transmitter  at  point  A  of  this  figure,  and  is  refracted  as  it  passes  through 
the  ionosphere,  where  it  meets  a  scattering  center  of  vertically-aligned  scatterers  at  point  B.  The 
scattered  rays  are  constrained  to  have  the  same  zenith  angle  as  the  incident  ray,  but  are  free  to 
have  any  azimuth  angle.  The  ray  reaching  the  receiver  is  thus  scattered  into  the  appropriate 
azimuth  in  which  it  continues  to  be  refracted,  eventually  turning  downward  so  it  eventually  leaves 
the  ionosphere  and  arrives  at  the  receiver  (point  D).  The  two  planes  of  propagation  defined  by  the 
transmitted  ray  and  the  vertical  scattering  center,  and  the  received  ray  and  scattering  center  meet 
at  the  line  BC  defined  by  the  scattering  center.  The  problem  can  be  simplified  by  considering  the 
two  propagation  planes  as  a  single  plane,  folded  about  the  line  BC.  This  single  plane  is  then 
flattened  so  that  the  whole  path  can  be  plotted  in  two  dimensions,  as  shown  in  Figure  4. 
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Figure  4  shows  the  unfolded  propagation  path  described  above,  along  with  the  vertical 
scattering  center.  The  preservation  of  elevation  angle  (or  its  complement  the  zenith  angle  ([))  in  the 
scattering  process  is  clear,  so  that  the  equation  of  the  scattered  path  is  the  same  as  that  of  a  ray 
transmitted  at  the  same  elevation  angle  and  refracted  (but  not  scattered)  by  the  ionosphere  so  it 
arrives  at  the  ground  at  a  point  separated  from  the  transmitter  by  a  horizontal  distance  which  is  the 
sum  of  the  distances  from  transmitter  to  scattering  center  and  scattering  center  to  receiver. 

The  transmitter  location  is  defined  by  its  (x,y)  coordinates  (x,,  y,) .  Likewise  the  receiver 
location  is  given  as  (x^,y^) .  For  a  scattering  center  located  at  (x^yg) ,  the  distances  cIq,  d-^  of 
Figure  4  are  found  from 

do= 

(6) 

dj  =  ( (x,  -  X,)  2  +  (y,  -  y,)  2) 
and  the  horizontal  distance  from  transmitter  to  receiver  is 

D  =  +  dj  C^) 


2.4  Modelled  Ionosphere 


The  ionosphere  is  modelled  as  a  single  F-layer  with  a  parabolic  electron  density  profile  is 
given  by 


N(h)  =No\l 
0 


V  Xm 


fork  =  hQ-y^to  hQ  +  y^  , 
elsewhere. 


(8) 


where  Nq  is  the  maximum  electron  density,  Hq  the  corresponding  height  of  maximum,  and  y^  the 
ionospheric  half  width,  shown  in  Figure  4.  The  refractive  index  p,  ignoring  collisions  and 
magnetic  field  effects  [12],  is  given  by 


(9) 


where  fp  is  the  plasma  frequency  given  by 


Ne^ 


(10) 


where  e,  m  are  the  electron  charge  and  mass,  and  is  the  permittivity  of  free  space.  From 
equations  (8),  (9),  and  (10),  the  refractive  index  p  can  be  expressed  in  terais  of  the  ionospheric 


parameters  by 

p2  ^  i_p2  +  j^2  ^ 


u  -^0  ^ 

where  9  ~  J  “ 


(11) 


for  heights  between  h^-y^  and  hQ  +  y^.  p  is  equal  to  1.0  elsewhere  (free  space).  The  critical 
frequency  /q  (which  is  the  plasma  frequency  at  the  height  of  maximum)  is  an  operator-entered 
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parameter,  along  with  the  height  of  maximum  Hq  and  ionospheric  half- width  y^. 


2.5  Propagation  Path  Calculation 

The  next  step  is  to  calculate  the  starting  zenith  angle  (t)Q  (see  Figure  4)  needed  for  a  ray  to 
propagate  to  a  receiving  array  at  ground  distance  D,  and  the  corresponding  distance  focusing 
aD/dtPo  for  that  distance.  The  laws  of  refraction  and  the  refractive  index  profile  given  by  (11)  are 
used  for  that  purpose. 

For  a  horizontally-stratified  ionosphere,  Snell's  law  states  that  for  any  height,  psin<t)  is 
constant,  where  ^  is  the  zenith  angle,  so  that 

psintj)  =  sin(t)Q  .  (12) 

The  horizontal  distance  along  the  ray  path  is  denoted  by  d.  Along  the  ray  path,  a  small  change  in 
horizontal  distance  bd  is  related  to  a  change  bh  in  vertical  distance  by 

bd  =  (tan(t))5/i.  (13) 


Substituting  (12)  in  (13),  we  get  sind)^ 

bd  =  ” - bh  . 

(sin(l)o)2 


Equation  (11)  may  be  substituted  in  (14)  and  the  result  integrated  from  the  transmitter  to  a 
horizontal  distance  d  less  than  or  equal  to  the  midpoint  distance  D/2,  to  obtain 


(^0  “ ^m)  tan(t)o  +  y  sin(|)olog 


X+  (cos(|)o)2-p2-i-x^ 
COS<t)Q  -I-  p 


At  the  midpoint  (Figure  4),  the  ray  is  horizontal  so  that  equation  (12)  becomes 

=  sin(l)o  ,  (16) 

where  p.  is  the  refractive  index  at  this  point.  Setting  d  =  D/2  ,  and  substituting  equations  (11) 
and  (16)  mto  (15),  the  ground  distance  D  is  found  as  a  function  of  the  starting  zenith  angle  (t)^  : 

Jl-a^  ym  L - 2,  + 

D  =  2ihQ-y„) — - —  +  log(^— ^)  where  a  =  cos<t)Q  .  (17) 

For  an  arbitrary  ground  distance  D  there  will  be  two  solutions  of  (17)  for  <j)Q  (high-angle  and  low- 
angle),  provided  that  any  exist  at  all. 

In  the  implemented  modelling  software,  a  stepwise  approach  is  taken  to  find  the  solutions 
for  <t)o  ,  and  at  the  same  time  to  find  the  distance  focusing  3D/3(pQ  .  The  search  is  start^  with  the 
ray  that  just  reaches  the  base  of  the  ionosphere  at  the  midpoint  distance  D/2.  The  zenith  angle  of 
this  ray  <t)Q  (start)  is  an  upper  bound  for  the  possible  zenith  angles  for  a  ray  propagating  to  a 
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ground  distance  of  Z).  It  is  given  by 


♦o  (18) 

A  check  is  made  to  see  if  the  ray  corresponding  to  this  zenith  angle  penetrates  the  ionosphere.  If  it 
does  then  all  lesser  zenith  angles  will  also  penetrate  the  ionosphere  and  there  will  be  no  solutions. 
Penetration  occurs  when  the  ray  can  reach  the  height  of  maximum  and  not  be  bent  horizontal,  i.e., 
when  (p  (/iq)  <  Tt/2.  From  (11)  and  (12),  penetration  occurs  when 

cos(t)Q>/o//  .  (19) 

Equation  (19)  is  used  to  check  for  the  existence  of  solutions. 

The  zenith  angle  is  then  incrementally  reduced  from  the  upper  limit  ^Qistart)  and  the 
‘search’  ground  distance  D'  found  for  each  step  from  equation  (17).  Reference  [13]  provides  a 
helpful  description  of  path  length  on  zenith  angle.  Initially  the  search  distances  D  are  greater 
than  the  actual  distance  D,  and  decrease  as  the  zenith  angle  is  reduced  (and  the  elevation  is 
raised).  When  the  search  distance  D'  first  becomes  less  than  the  actual  ground  distance  D,  that 
value  and  its  corresponding  zenith  angle  and  the  value  for  the  previous  step  are  interpolated,  to 
find  the  zenith  angle  where  D'  =  D.  The  interpolation  takes  the  form  of  an  iterative  linear 
interpolation,  with  the  change  in  zenith  angle  decreasing  each  iteration,  until  it  reaches  a  limit  set 
by  the  program  (10'^  radians).  The  zenith  angle  thus  obtained  is  the  low-angle  solution  (low) . 
The  distance  focusing  parameter  9Z)/9(Pq  (low)  is  also  found  by  this  interpolation  proc^ure. 

The  stepping  downward  in  zenith  angle  (i.e.,  increasing  in  elevation)  is  then  continued.  The 
search  distances  D'  continue  to  decrease  for  a  number  of  steps,  to  a  minimum  value  (the  skip 
distance)  after  which  they  begin  to  increase.  The  point  at  which  D'  again  becomes  greater  than  Z) 
is  noted.  That  point  and  the  previous  point  are  interpolated  as  before  to  find  the  high-angle  zenith 
angle  (|)q  (high)  and  focusing  parameter  BD/d^  (high) . 

It  should  be  noted  that  when  the  transmitted  frequency  is  well  below  the  maximum  useable 
frequency  (MUF)  for  the  distance  D,  the  zenith  angle  for  the  high-angle  ray  will  be  less  than  one 
step  above  0°.  The  step-search  procedure  will  fail  in  such  cases,  as  a  step  decrease  in  the  zenith 
angle  will  reduce  it  to  a  value  less  than  zero,  which  is  physically  impossible.  These  cases  can  be 
safely  ignored,  as  the  focusing  parameter  dD/d(pQ  becomes  extremely  large  for  the  high-angle 
ray  at  frequencies  well  below  the  MUF  and,  as  will  be  seen  in  the  next  section  (equation  (28)), 
this  results  in  a  very  low  signal  strength.  In  such  cases  the  high-angle  ray  is  not  included. 


2.6  Signal  Strength  Calculation 

2.6.1  Mean  Received  Power  for  a  Ray  Path 

For  each  valid  path,  the  simulation  program  computes  the  mean  signal  power,  and  a 
sequence  of  signal-sample  contributions  at  the  receiving  array.  This  calculation  depends  on  the 
geometry,  ionospheric  focusing,  and  scattering  efficiency. 
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center 


Figure  5.  Transmitter-to- 
scattering  center  geometry 
showing  subtended  angles. 


A  single  vertical  scattering  center  may  be  considered  to  have  a  cross  section  per  unit  length 
or  ‘scattering  width’  w,  which  can  be  treated  as  a  horizontal  distance.  As  shown  in  Figure  5,  it 
subtends  an  angle  a  at  the  transmitter,  given  by 


(20) 


The  transmitted  power  radiated  into  this  angle  will  be  reradiated  by  the  scattering  center.  For  an 
incremental  change  in  zenith  angle  d(^Q  at  the  transmitter,  the  corresponding  solid  angle  which 
includes  the  scattering  center  is 

wsin(t)Q 

dO.  =  asin(t)od<l)o  =  — — d^Q.  (21) 

“0 


The  power  radiated  into  dCl  is 


dP, 


Pfi^dQ. 

47C 


(22) 


where  F,,  G,  are  the  transmitter  radiated  power  and  antenna  gain  respectively.  Combining  (21) 
and  (22),  the  following  is  obtained: 


dP 


S 


PjG^ivsin<l)Q 

A’TzcIq 


(23) 


Note  that  the  power  dP^  will  be  reradiated  by  the  scattering  center  equally  into  all  azimuths, 
with  no  change  in  the  zenith  angle.  The  scattering-center-to-receiver  geometry  is  illustrated  in 
Figure  6.  At  the  receiving  array,  the  change  in  zenith  angle  d^Q  corresponds  to  a  horizontal 


distance 

dD  =  |8D/a(t)o|d(|)o  . 


(24) 
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Figure  6.  Scattering  geometry  showing  region  of  coverage  of  scattered  signal  energy. 

The  corresponding  distance  perpendicular  to  the  ray  direction  is 

d5  =  dDcos^Q  ,  (25) 

since  the  symmetrical  nature  of  the  ionospheric  refraction  along  the  path  implies  equal 
transmitting  and  receiving  zenith  angles  <j)Q.  Tiie  scattered  energy  passes  through  a  ring  of  width 
db  and  diameter  d^  centered  on  the  point  of  intersection  of  the  scattering  center  with  the  ground. 
The  area  of  this  ring,  normal  to  the  ray  path,  is 

dA  =  2%d^db  =  2jtdiCOS<l)Q|8D/3<l)Q|d(|)Q  . 

From  (23)  and  (26),  the  power  density  at  the  receiving  array  is  found  to  be 

_dP^  ^  P,G,w tan(l>o 

~dA  ~  87t2dodi|aZ)/a(t>o|  ' 

The  received  power,  for  an  array  element  antenna  of  gain  G^,  is  then 

S^G^X'^  PfGp^X^wtm^Q 
47U  32Tt^dod^\dD/d<^Q\  • 


(26) 

(27) 

(28) 
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It  remains  to  specify  the  effective  scattering  width  w  as  a  function  of  height  and  time,  before 
the  rms  signal  strength  can  be  obtained. 


2.6.2  Effective  Scattering  Width 


The  physical  processes  involved  in  the  actual  radio  wave  scattering  determine  the 
dependence  of  effective  scattering  width  w  upon  the  excess  electron  density  &N  and  size  of  the 
scattering  irregularities.  A  radio  wave  incident  on  a  vertical  field-aligned  irregularity,  of 
horizontal  area  A  corresponding  to  small  transverse  dimensions  relative  to  a  wavelength,  will  by 
virtue  of  its  electric  field  induce  an  excess  oscillating  vertical  current  in  the  irregularity 
proportional  to  the  excess  electrons  per  unit  length,  given  by  A6iV.  The  scattered  wave  can  be 
considered  to  be  radiated  by  this  excess  current.  As  such,  it  will  have  an  associated  magnetic  field 
B  proportional  to  the  excess  current,  and  therefore  to  AbN,  and  a  corresponding  electric  field 
proportional  to  the  rate  of  change  of  current  which^is  also  proportional  to  AdN.  The  power 
density  of  the  scattered  wave  is  proportional  to  ExB  and  therefore  (AbN)^  by  the  preceding 
arguments.  In  the  present  model,  the  density  and  sizes  of  the  vertical  scatterers  in  a  particular 
region  are  assumed  to  be  constant  over  the  region,  and  the  scattering  centers,  likewise  constant,  so 
that  each  center  can  be  said  to  have  an  equal  number  of  scatterers  similar  in  their  size  distribution. 
Under  these  assumptions,  the  scattered  wave  from  a  collection  of  scatterers  comprising  a 
scattering  center  will  have  a  mean  power  density  proportional  to  the  weighted  mean  square  excess 
electron  density  of  those  scatterers.  The  effective  scattering  width  for  the  scattering  center  must 
also  be  proportional  to  the  mean  square  excess  electron  density. 

f  . 

In  order  to  relate  the  scattering  width  w  to  the  height  of  ^  scattenng,  two  options  have 
been  included.  In  the  first  option,  the  rms  excess  electron  density  dN  is  taken  to  be  proportional  to 
the  electron  density  N.  (This  assumption  has  been  used  by  other  researchers,  as  being  the  simplest 
assumption  in  the  absence  of  experimental  data  [14].)  Under  this  assumption,  the  scattering  width 
w  is  proportional  to  N^.  From  (8),  w  can  be  written  as  a  function  of  height,  by 


W  =  Wq 


(29) 


where  Wq  is  the  effective  scattering  width  at  the  hei^t  of  maximum  electron  density  Hq.  In  the 
second  option,  8N  is  assumed  to  be  proportional  to  Jn,  so  that  w  is  proportional  to  N.  This  gives 
a  somewhat  less  height-dependent  scattering  width,  of  the  form 


W  =  Wn 


1- 


The  proportionality  constant  Wq  is  an  operator-specified  parameter. 


(30) 


It  remains  to  find  the  height  h,  of  the  scattering,  before  w  can  be  obtained.  Under  the 
scattering-center  assumption,  the  scattering  effectively  takes  place  at  a  single  height, 
corresponding  to  the  location  of  the  scattering  center.  (This  is  a  simplification  of  what  actually 
happens,  since  in  the  real  world,  the  individual  scatterers  will  be  some  distance  apart,  and  the  ray 
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For  an  identified  ray  path  for  which  the  initial  zenith  angle  <1)q  has  been  obtained,  equations  (32), 
(31),  and  (33),  in  that  order,  are  used  to  obtain  the  scattering  height.  Since  the  ray  pth  is 
symmetric  about  point  of  maximum  height  in  the  unfolded  plane  of  propagation,  and  the  distance 
for  which  these  formulas  apply  must  be  less  than  D/2,  the  lesser  of  do 
taken  as  d  in  (32).  The  effective  scattering  width  is  then  calculated  firom  (29)  or  (30)  according  to 
the  user’s  requirement. 

The  simulation  program  then  calculates  the  mean  received  power  for  the  ray  path  from  (28). 
Included  in  this  calculation  are  the  antenna  gains  G,  and  which  can  be  specified  by  the  user  as 
a  function  of  zenith  and  elevation  angle,  in  the  simulation  software. 


2.6.3  Generating  a  Sequence  of  Signal  Samples 

For  each  scattering  path,  a  sequence  of  signal-sample  contributions  will  have  the  mean 
signal  power  calculated  above,  from  which  an  rms  amplitude  can  be  obtained.  However,  the 
individual  sample  contributions  will  be  affected  in  their  phase  and  amplitude  by  sample-to-sample 
changes  in  the  individual  scatterers  making  up  the  scattering  center.  These  changes  are  simulated 
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by  employing  bandwidth-limited  white  noise  where  the  bandwidth  is  selected  to  be  similar  to  the 
Doppler  spreading  observed  in  practice  over  high-latitude  propagation  paths  (e.g.,  10  or  20  Hz), 
accordingtothefollowingprocedureprovidedby  Jorgenson[15]:  The  simulation  program  upon 
’  its  initiation  generates  the  appropriate  digital  filter  coefficients  for  a  fourth-order  Gaussian 

spectral  shaping  filter,  according  to  the  user-specified  Doppler  bandwidth  and  sampling  rate. 
Later,  for  each  path,  it  generates  a  sequence  of  completely  random  complex  samples  of  unit  rms 
^  amplitude,  and  passes  this  sequence  through  the  filter  to  obtain  the  required  bandwidth-limited 

noise  samples.  These  complex  samples  are  multiplied  by  the  rms  amplitude,  to  generate  a 
sequence  of  sample  contributions  for  the  path  in  question,  as  seen  at  the  receiver  reference  point. 

The  sample  contributions  are  then  adjusted  in  their  phase,  according  to  the  specified  x,  y 
position  of  each  array  element  and  the  incoming  zenith  and  azimuth  angles,  to  obtain  a  sequence 
of  sample  contributions  for  each  array  element.  The  signal-sample  contributions  are  accumulated 
over  all  paths,  on  a  sample-by-sample  basis,  for  each  of  the  array  elements,  to  generate  a  set  of 
signal  samples  seen  by  the  array,  due  to  scatter  propagation. 


2.7  Sporadic-E  Propagation 


Sporadic  E  (Eg)  propagation  is  simulated  by  assuming  a  horizontal  mirror-like  reflection  at  a 
height  hg,  midway  between  the  transmit  and  receive  locations  along  the  great-circle  meridian. 
The  value  used  for  hg  is  100  km,  which  is  typical  of  Eg  heights.  A  user-specified  reflection 
coefficient  R,  representing  the  fraction  of  signal  power  reflected,  is  used  to  specify  the  reflection 
efficiency. 


The  total  path  length  L, 
geometry  to  be 

L  = 


for  the  modelled  Eg  propagation  is  determined  by  the  flat-earth 
ax,-x,)^*{y-y,)U4hl)''^  (34) 


and  the  corresponding  received  power  is 


(35) 


The  corresponding  zenith  and  azimuth  angles  <]),  0  are 

2/z 

<t)  =  acos(-^)  and  0  =  atan  ( (x,-x^),  •  (36) 

The  received  Eg  signal  is  assumed  to  be  constant  in  amplitude  and  phase  between  samples.The 
amplitude  is  obtained  from  the  receive  power,  which  is  calculated  from  equations  (34),  (36),  and 
(35).  The  phase  is  chosen  randomly  for  the  first  sample  and  maintained  constant  thereafter.  The  Eg 
contribution  at  each  element  of  the  array  is  then  computed  by  adjusting  the  phase,  according  to  the 
element  x,  y  position  and  the  signal  direction.  This  contribution  which  is  constant  over  all 
samples,  is  added  to  the  scattered  signal  seen  by  each  element,  to  get  the  total  received  signal  for 
that  element. 


15 


3.0  DIRECTION-FINDING  SYSTEM  MODEL 


The  receiving  array  employed  by  the  simulation  package  can  have  a  variable  number  of 
element  antennas,  with  arbitrary  x,  y  positions,  measured  relative  to  the  receiving  reference  point 
X  ,  y  .  At  the  time  of  writing,  two  options  exist  for  the  receiving  and  transmitting  element  antenna 
patterns:  an  isotropic  radiation  pattern,  and  a  cos(elevation)  pattern  similar  to  that  experienced  by 
a  short  whip  antenna  over  perfectly  conducting  ground. 

Two  techniques  have  been  employed  in  the  simulation  package,  to  perform  the  direction¬ 
finding:  a  plane-wave  fitting  technique,  and  the  MUSIC  algorithm. 


3.1  Plane- Wave  Fitting  Procedure 


An  iterative  search  procedure  is  invoked  to  find  the  best  plane-wave  fit  to  the  total  received 
signal  at  the  array  elements.  The  signals  seen  by  the  N  element  antennas  can  be  viewed  as  the 
complex  components  of  a  iV-dimensional  ‘signals’  vector.  A  hypothetical  plane  wave  would 
likewise  have  its  corresponding  N-dimensional  array  response  vector  which  could  be  calculated 
from  the  plane-wave  direction  and  amplitude.  The  best  plane-wave  fit  has  the  minimum  error 
power,  i.  e.,  the  length  of  the  difference  vector  between  the  plane  wave  response  vector  and  the 
total  received  signal  vector  seen  at  the  array  elements  is  a  minimum  over  all  plane  wave  directions 
and  amplitudes. 

It  can  be  shown  that  the  array  signals  vector  calculated  for  a  unit-amplitude  plane  wave  in  a 
given  direction  will  have  a  correlation  C  with  the  received  signals  vector  which  is  a  maximum 
when  the  direction  is  the  best-fit  direction.  The  correlation  is  given  by 

c  =  i  v'-  (37) 

n  =  1 

where  is  the  complex  received  signal  sample  at  the  nth  antenna,  and  6^  the  phase  of  the  plane 
wave  at  the  nth  antenna,  which  is  given  in  terms  of  the  plane  wave  bearing  a  and  zenith  angle  p, 
by 

5  = — ^  (x„sina-i-y„cosa)  cosP  ,  (38) 

”  K 

where  x  ,y  are  the  antenna  coordinates.  The  corresponding  best-fit  amplitude  is  given  by 

A  =  C/N  (39) 

The  iterative  search  procedure  described  below  is  used  to  find  the  best-fit  plane-wave 

direction  (i.e.,  bearing  a  and  zenith  angle  p): 

1.  First,  the  number  of  sidelobes  that  could  exist  in  the  pattern  of  the  array  when  steered  to  a 
direction  are  estimated,  by  considering  the  total  change  in  phase  that  can  occur  between  the 
furthermost  elements  of  the  array,  in  going  from  signals  pa^lel  and  antiparallel  to  a  particular 
axis.  These  numbers,  for  the  x  and  y  axes  correspond  to  the  number  of  local  maxima  that  can 
exist  in  the  correlation  function  C  over  all  directions. 
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2.  Next  a  suitable  number  of  steps  in  direction  cosine  are  selected,  to  ensure  that  each  sidelobe 
or  local  maximum  will  be  seen.  In  the  model,  the  number  of  steps  riy  were  set  to  5  times 
the  number  of  possible  sidelobes,  in  the  x  and  y  directions.  The  direction  cosines,  in  the  x  and 
y  directions,  are 

=  -sinacosP  and  =  -cosacosp ,  (40) 

respectively.  The  direction  cosines  \  and  A^  are  stepped  from  0  to  1  in  steps  of  1/n^  and 
1  /n  ,  and  the  correlation  function  C  evaluated  for  each  combination  of  A^  and  A^ . 

3.  The  J  largest  values  of  C  are  selected,  along  with  their  direction  cosines,  as  possible  choices 
for  an  approximate  maximum. 

4.  For  each  of  these  points,  an  iterative  search  procedure  is  followed,  to  find  the  local 
maximum.  This  consists  of  halving  the  step  size,  and  computing  C  for  A^  and  A^  a  step  away 
from  the  present  point,  in  both  x  and  y  directions.  These  new  values  are  compared  with  the 
previous  C-value,  and  the  point  corresponding  to  the  largest  value  chosen.  The  step  size  is 
halved  again,  and  the  procedure  repeated.  The  process  is  continued  until  the  step  size  falls 
below  a  predetermined  small  value  (.01).  The  final  point  and  its  C-value  are  taken  as  the  local 
maximum. 

5.  The  largest  of  the  8  local  maxima  thus  found  is  taken  as  the  maximum.  The  corresponding 
best-fit  direction  is  found  from  the  corresponding  direction  cosines  A^  and  A^  by 

-A  A 

a  =  atan  f— A  —A  )  ,  B  =  acos  (—:——)  or  B  =  ( - )  (41) 

-  duun,  LXy)  ,  H  ^sina  cosa 

The  corresponding  complex  amplitude  of  the  best-fit  plane  wave  is  found  from  equation  (39). 

The  difference  between  the  total  signal  at  each  of  the  array  elements  and  the  corresponding 
best-fit  plane  wave  signal  is  estimated.  The  power  of  this  difference  is  summed  across  the 
antennas  and  compared  to  the  total  received  power  also  summed  across  the  antennas.  The 
resultant  error-power  fraction  is  reported  by  the  plane-wave-fitting  program,  as  a  quality  factor 
representative  of  the  goodness  of  fit. 

The  plane-wave-fitting  procedure  requires  only  a  single  sample,  and  provides  a  single 
representative  arrival  direction  as  well  as  a  plane-wave  goodness-of-fit  value. 


3.2  MUSIC  Algorithm 


This  section  provides  a  brief  theoretical  discussion  of  the  MUSIC  algorithm  [11],  and  a 
description  of  how  it  is  implemented.  r  -i 


(0 


The  signals  received  by  the  N-element  array  can  be  written  as  a  vector  x{t)  = 


%(0 


^  ^// 

where  (r)  denotes  time  dependence.  The  covariance  matrix  of  these  signals  is  given  by  R  =  xx  , 
where  the  averaging  is  denoted  by  the  overhead  bar  ,  and  a  conjugate  transpose,  by  the 
superscript  In  order  to  compute  R,  the  averaging  is  taken  over  a  sequence  of  time-samples 
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generated  by  the  simulation.  The  normalized  eigenvectors  of  R  are  denoted  by  e^, 
arranged  according  to  the  size  of  their  corresponding  real  eigenvalues  Xj, in  increasing 
order. 

For  5  independent  signals  incident  on  the  array  ( 5  <  ),  the  eigenvectors  _  5  +  j , . . % 
span  the  signal  space  generated  by  the  signals,  while  the  remaining  (orthonormal)  eigenvectors 
span  the  noise  space. 


A  unit-amplitude  plane- wave  signal  vector  r  (a,  p)  corresponding  to  a  particular  direction, 

i.e.  azimuth  a  and  zenith  angle  P,  has  the  form 


r(a,  p) 


_1_ 

JN 


iS  (aP) 
e 

^»5^(a,P) 


(42) 


where  the  5^  (a,  P)  are  given  by  equation  (38).  For  a  direction  (a,  p)  to  correspond  to  a  signal 
in  the  signal  space,  its  corresponding  signal  vector  r  (a,  P)  must  be  perpendicular  to  the  noise 
space.  That  is,  its  projection  length  onto  the  noise  space,  given  by 


n(a,  P) 


(o,p)ir 


1/2 


(43) 


is  zero.  In  practice,  a  minimum  will  do.  The  MUSIC  algorithm  consists  of  making  a  guess  at  the 
number  of  ‘signal’  directions  S',  then  using  the  S'  eigenvectors  corresponding  to  the  S'  lowest 
eigenvalues  to  define  the  noise  space  (i.e.,  setting  S  =  S').  Equations  (42)  and  (43)  are  then  used 
along  with  a  search  procedure,  to  finding  the  S'  directions  (a,  P)  corresponding  to  the  S  lowest 
local  minima  of  11  (a,  P) . 

The  implemented  procedure  for  determining  signal  directions  according  to  the  MUSIC 
algorithm  is  as  follows: 

1.  Find  the  covariance  matrix  R  for  the  sequence  of  signal  samples,  and  from  it  calculate  the 
corresponding  normalized  eigenvectors  Cj, ...,  C/v  arranged  according  to  increasing 
eigenvalue.  (The  sequence  length  M^N  for  R  to  be  nonsingular.) 

2.  Specify  the  desired  number  of  directions  (S' )  to  search  for. 

3.  Searching  over  all  directions  (a,  P) ,  compute  11  (a,  P)  from  (38),  (42),  and  (43),  looking  for 
the  S'  lowest  local  minima.  The  search  technique  is  conducted  over  the  direction  cosines,  in  a 
manner  similar  to  that  described  in  section  3.1  for  the  plane- wave  fitting  technique.  The 
estimated  signal  directions  are  those  of  the  S'  lowest  minima.  The  corresponding  projection 
lengths  n  are  an  estimate  of  how  well  each  of  the  directions  fit  the  signal  space. 

The  implemented  MUSIC  algorithm  requires  a  sequence  of  signal  samples  from  each  of  the 
array  elements,  and  an  operator  request  as  to  the  number  of  signal  directions  to  look  for.  It 
provides  the  requested  number  of  signal  directions  along  with  an  estimate  of  how  well  each 
direction  fits  into  the  signal  space  represented  by  the  sequence  of  signal  samples. 
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4.0  TESTS  OF  THE  FIELD-ALIGNED  SCATTERING  MODEL 

WITH  PLANE- WAVE  FITTING 

4.1  Operational  Parameters 

A  number  of  simulations  were  performed  to  examine  the  implications  of  the  field-aligned 
scattering  model,  and  to  compare  the  results  with  those  of  the  previously  acquired  Kestrel  data 
[9,10].  These  simulations  used  a  receiving  array  geometry  and  a  transmitting  location  similar  to 
that  of  the  Kestrel  measurements  on  the  Thule  transmissions.  The  receiving  array  consisted  of 
seven  elements  at  specified  locations  about  a  circle  of  25  m  radius.  The  transmitting  antenna  was 
assumed  to  be  700  km  distant,  and  the  operating  frequency  was  set  to  8.0  MHz.  The  transmitting 
and  receiving  antenna  patterns  were  assumed  to  be  those  of  a  short  vertical  whip  on  a  perfectly 
conducting  ground  plane.  Although  these  patterns  do  not  match  those  of  the  Kestrel  experiment, 
they  do  possess  the  overhead  null  of  the  actual  antenna  patterns  and  are  sufficient  to  demonstrate 
the  model.  The  assumed  transmitted  power  was  1  kW,  and  the  receiver  threshold,  -115  dBm. 

The  ionosphere  was  assumed  to  have  a  F-region  height  of  maximum  of  350  km,  and  a 
halfwidth  of  100  km.  Various  values  were  attempted  for  the  critical  frequency  /q,  from  5  to  10 
MHz.  Scattering  region  or  patch  sizes  were  200  by  200  km.  The  scattering-center  density  inside 
these  patches  was  taken  as  constant;  values  for  different  tests  ranged  from  .01  to  0.1  scattering 
centers/km^  corresponding  to  400  to  4000  scattering  centers  per  patch,  respectively.  The  effective 
scattering  width  Wq  at  the  height  of  maximum  Hq  was  arbitrarily  selected  to  be  10  m. 

The  plane-wave  fitting  DF  technique  was  employed  for  direction  estimation.  This  technique 
was  similar  in  most  respects  to  the  technique  used  for  the  Kestrel  data. 


4.2  Simulation  Results 
4.2.1  Single  Scattering  Patch 

A  series  of  simulations  were  run,  using  the  simple  case  of  a  single  scattering  patch  with  a 
horizontal  motion  perpendicular  to  the  direct  transmitter-receiver  path  (referred  herein  as  the 
great-circle  path,  despite  the  flat-earth  geometry) ,  crossing  the  great-circle  path  at  the  midpoint. 
This  geometry  is  illustrated  in  Figure  7.  A  2-km  shift  of  the  scattering  patch  was  used  between 
estimates,  which  corresponds  to  a  4-second  time  for  a  drift  velocity  of  500  m/sec  which  is  typical 
of  the  polar  cap. 

QTx 

patch  of 
scatterers 


Figure  7.  Plan  view  of 
single-patch  scattering 
geometry. 
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In  the  first  test,  a  critical  frequency  /q  of  7.0  MHz  was  used.  This  corresponds  to  a  situation 
where  the  corresponding  ionospheric  electron  density  is  a  little  higher  than  the  minimum  needed 
to  support  propagation  of  the  8  MHz  transmitted  signal  along  the  700  km  great-circle  transmitter- 
receiver  path.  The  density  of  scattering  centers  was  assumed  to  be  .02/km^.  The  stronger  of  the 
two  postulated  scattering-width  height  dependencies  (equation  (29)  was  used,  corresponding  to 
the  excess  electron  density  of  the  scatterers  being  proportional  to  the  ambient  electron  density. 

Figures  8a,  b,  c,  and  d  show  the  results  obtained.  Figures  8a  and  b  show  the  estimated 
plane-wave  bearings  and  elevations  as  a  function  of  estimate  number  or  time.  (A  0°  bearing 
corresponds  to  the  great-circle  direction  from  the  transmitter.)  The  bearings  and  elevations  appear 
to  follow  the  location  of  the  patch  of  scatterers,  with  apparently  random  point-to-point  variations. 

The  point-to-point  bearing  variations  cover  a  range  of  values  roughly  equal  to  the  bearing 
angle  subtended  by  the  patch  at  the  receiving  array,  at  any  one  time.  The  elevation  variations,  on 
the  other  hand,  are  substantially  less  than  the  elevation  angle  subtended.  For  instance,  at  the 
midpoint  of  the  simulation  when  the  patch  is  halfway  between  the  transmitter  and  receiver,  it 
subtends  a  bearing  angle  of  33®  centered  on  0°,  and  at  the  height  of  maximum  electron  density, 
an  elevation  angle  of  17°  centered  on  45°.  An  inspection  of  Figures  8a  and  b  for  this  time 
indicate  a  similar  bearing  spread  for  the  direction  estimates,  but  a  much  smaller  elevation  angle 
spread.  The  elevation  angles  with  a  few  exceptions  lie  within  a  5°  range  centered  on  43°.  This  is 
likely  due  to  the  fact  that  the  strongest  contributions  to  the  signal  come  from  those  scattering 
centers  which  intersect  the  generally  better-focused  low-angle  ray  nearest  the  height  of  maximum 
electron  density.  This  is  a  result  of  the  larger  scattering  cross  section  as  well  as  the  increased 
focusing  that  occurs  for  low-angle  rays  which  approach  those  heights.  These  scattering  centers  are 
located  in  the  patch  along  the  bisector  of  the  transmit-receive  great-circle  path,  where  the 
elevation  angles  for  scatter  propagation  are  expected  to  lie  in  a  narrow  range  near  45° . 

Figure  8c  shows  the  estimated  signal  directions  obtained  in  the  run,  in  an  elevation-versus- 
bearing  plot.  The  directions  are  scattered  about  a  well-defined  arc  centered  on  the  great-circle  (0°) 
bearing.  This  behaviour  is  expected,  since  ray  paths  closer  to  the  great-circle  path  will  by  virtue  of 
their  shorter  distances  correspond  to  higher  elevation  angles  which  penetrate  higher  into  the 
ionosphere. 

Figure  8d  is  a  time  history  of  the  number  of  scattering  points,  or  propagation  paths 
contributing  to  the  signal  seen  by  the  array.  With  the  specified  density  of  scattering  centers  (.02/ 
km^)  and  the  scattering  patch  size  (200  km  square),  there  are  800  scattering  centers,  each  possibly 
contributing  a  high-angle  and  a  low-angle  ray  to  the  signal.  For  most  of  the  time,  the  number  of 
scattering  points  in  Figure  8d  is  slightly  less  than  800.  This  implies  that  only  the  stronger  low- 
angle  rays  contributed  to  the  signal.  The  high-angle  rays  for  most  scatterers  had  near-zero  zenith 
angles  for  which  the  defocusing  was  very  high,  and  so  were  rejected  by  the  simulation  program. 
At  the  midpoint  of  the  run,  the  number  of  contributing  paths  rises  slightly  above  800.  The  patch  is 
midway  between  the  transmitter  and  receiver  at  this  time,  so  that  the  7.0  MHz  transmitter 
frequency  was  closer  to  the  MUF  (maximum  useable  frequency)  for  the  scatter  propagation  paths. 
The  high-angle  ray  under  these  circumstances  would  not  lie  as  close  to  the  zenith,  so  that  its 
contribution  would  be  included. 
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The  resultant  signals  were  observed  to  fit  a  single  plane  wave  remarkably  well.  Of  the  500 
estimates  represented  in  Figures  8a,  b,  c,  and  d,  398  had  an  error  power  fraction  (Section  3.1)  of 
less  than  or  equal  to  0.1.  In  other  words,  for  398  out  of  500  estimates,  90  percent  or  more  of  the 
received  signal  power  could  be  attributed  to  a  single  plane  wave. 


Figure  8a.  Estimated 
plane-wave  bearing, 
as  a  function  of 
estimate  number,  for 
the  single-patch 
geometry  of  Figure  7 
and  a  critical  frequency 
of?  MHz. 


200  300 

TIME  (estlmoie  no.) 


Figure  8b.  Estimated 
plane-wave  elevation, 
as  a  function  of 
estimate  number. 


Figure  8c.  Estimated 
plane-wave  directions, 
plotted  in  terms  of 
bearing  and  elevation. 
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4.2.1.1  Dependence  on  Ionospheric  Electron  Density 

The  ionospheric  electron  density  is  specified  in  terms  of  the  critical  frequency  /q  ,  and 
determines  the  maximum  firequency  (MUF)  that  will  propagate  between  any  two  points  a 
specified  ground  distance  apart.  In  testing  the  model,  /q  was  varied  to  obtain  great-circle  MUFs 
that  were  a  little  greater  than,  nearly  equal  to,  and  less  than,  the  transmitting  frequency  /.  These 
values  of /o  were  7.0,  6.5,  and  6.0  MHz  respectively.  The  geometry  was  that  of  Figure  7,  and  the 
remaining  ionospheric  and  scattering  parameters,  the  same  as  previously  used.  The  results  are 
shown  in  Figures  9  and  10. 


TIME  (estimate  no.) 


fo  =  7.0  MHz 
(transmitted 
frequency  less 
than  great-circle 
MUF) 


fo  =  6.5  MHz 
(transmitted 
frequency  near 
great-circle 
MUF) 


fo  =  6.0  MHz 

(transmitted 

frequency 

greater  than 

great-circle 

MUF) 


Figure  9.  Number  of  contributing  paths,  as  a  function  of  estimate  number,  for  the 
single-patch  geometry  of  Figure  7  and  various  critical  frequencies  fo. 
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BEARING  (dagraat)  _  BEARING  (dagraaa)  _  ___  BEARING  j;dagraaa) 


fo  =  7.0  MHz  (transmitted  frequency  less  than  great-circle  MUF) 


fo  =  6.0  MHz  (transmitted  frequency  greater  than  great-circle  MUF) 


Figure  10.  Estimated  plane-wave  bearing  and  elevation,  as  a  function  of  estimate  number, 
for  the  single-patch  geometry  of  Figure  7  and  various  critical  frequencies  fo. 


Figure  9  shows  the  number  of  scatter  points  or  paths  contributing  to  the  signal,  as  the  patch 
moves  from  right  to  left  of  the  transmitter-receiver  great-circle  path,  for  the  three  critical 
frequencies 

The  behaviour  for  U  =  7.0  MHz  has  already  been  discussed.  In  this  case,  the  MUF  for  the 
great-circle  path  is  above  the  transmitted  frequency  so  propagation  over  this  path  is  supported. 

At  fo  =  6.5  MHz,  the  MUF  for  the  great-circle  path  drops  close  to  the  8.0  MHz  transnutted 
ftequency,  so  that  the  number  of  contributing  paths  begins  to  depend  strongly  on  patch  position. 
At  the  start  and  end  of  the  simulation  the  patch  is  well  removed  from  the  great  circle  so  that  the 
scattering  path  lengths  are  much  longer  and  the  MUFs  well  above  8.0  MHz.  The  number  of 
contributing  paths  is  near  800,  corresponding  to  one  path  per  scattering  center  (the  low-angle 
path).  The  high-angle  paths  are  dropped  as  a  result  of  their  being  very  close  to  the  zenith  and  thus 
extremely  defocused.  As  the  patch  approaches  the  great  circle,  the  MUFs  for  the  scattering  paths 
drop  toward  the  transmitted  frequency.  The  zenith  angles  for  the  high-angle  paths  increase,  and 
they  begin  to  be  included  in  the  signal  estimation.  The  number  of  paths  seen  near  the  midpoint  of 
the  run  is  approximately  16(X)  as  a  result  of  both  the  low  and  high- angle  paths  of  all  scattering 
centers  in  the  patch  contributing  to  the  signal. 

At  fo  =  6.0  MHz,  the  great-circle-path  MUF  has  dropped  below  the  8.0  MHz  transmitted 
frequency  and  the  behaviour  is  radically  different.  Propagation  becomes  impossible  for  those 
scattering  centers  whose  horizontal  path  length  corresponds  to  a  MUF  below  8.0  MHz.  Thus  the 
number  of  contributing  paths  drops  as  the  patch  approaches  the  great  circle,  from  1600  to  zero. 

The  corresponding  bearing  and  elevation  behaviour  is  shown  in  Figure  10. 

As  the  critical  frequency  is  dropped  from  7.0  to  6.5  MHz,  there  is  little  change  in  the 
estimated  bearings.  A  noticeable  increase  is  seen  in  the  estimated  elevation  angles,  especially 
when  the  patch  approaches  the  great  circle.  This  is  due  to  the  added  contribution  of  the  high-angle 
rays  to  the  signal  as  well  as  the  increased  elevation  angles  for  the  low-angle  rays  which  are  able  to 
penetrate  further  into  the  less  dense  ionosphere. 

When  the  critical  frequency  drops  to  6.0  MHz,  below  that  needed  to  support  great-circle- 
path  propagation,  there  is  a  significant  change  in  the  estimated  bearings.  There  are  no  bearing 
estimates  for  the  times  where  the  patch  is  sufficiently  close  to  the  great  circle  that  none  of  the 
scattering  distances  support  propagation.  For  slightly  farther  patch  positions,  the  bearings  reflect 
the  positions  of  the  farthest  scattering  centers  in  the  patch,  whose  horizontal  propagation  paths  are 
long  enough  to  support  propagation.  When  the  patch  is  sufficiently  distant  from  the  great  circle, 
all  scattering  centers  support  propagation  and  the  bearings  are  centered  roughly  toward  the 
midpoint  of  the  patch. 

A  different  aspect  of  the  influence  of  the  ionosphere  on  scattering  propagation  is  seen  when 
a  critical  frequency  well  above  the  great-circle-path  MUF  is  used.  Figures  11a,  b,  and  c  show  the 
time  history  for  a  critical  frequency  of  10.0  MHz.  Other  model  parameters  are  identical  to  those 
used  previously. 
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Figure  11a  shows  the  number  of  scattering  paths  contributing  to  the  signal  as  a  function  of 
time,  as  the  patch  moves  from  right  to  left  of  the  great  circle.  The  number  (~3(X))  is  well  below  the 
number  of  scattering  centers  (~800)  in  the  patch,  and  remains  so  over  the  path  of  travel  of  the 
patch.  This  happens  because  the  portion  of  the  propagation  path  lying  inside  the  ionosphere  is 
quite  short  when  the  transmitted  frequency  is  well  below  the  MUF.  A  scattering  center  contributes 
to  the  propagation  only  when  its  propagation  paths  intersect  it  at  a  point  inside  the  ionosphere. 
The  region  where  the  potential  paths  lie  inside  the  ionosphere  at  the  point  of  scattering  is  in  this 
case  a  narrow  strip  which  crosses  the  great-circle  path  perpendicularly  at  its  midpoint.  Only  the 
scattering  centers  inside  this  region  can  contribute  to  the  signal.  The  number  of  contributing  paths 
at  the  midpoint  (~270)  suggest  that  this  region  is  roughly  270/800  x  200  km  =  67.5  km  wide  at  the 
point  where  it  crosses  the  great  circle. 


Figure  11a.  Number 
of  propagation  paths 
contributing  to  the 
signal,  as  a  function 
of  estimate  number, 
for  the  single-patch 
case  of  Figure  7  and 
a  critical  frequency 
of  10  MHz. 


Figure  11b.  Estimated 
plane-wave  bearing, 
as  a  function  of 
estimate  number,  for 
a  critical  frequency  of 
10  MHz. 


Figure  11c.  Estimated 
plane-wave  elevation, 
as  a  function  of 
estimate  number,  for 
a  critical  frequency  of 
10  MHz 


Figures  11b  and  c  show  the  bearing  and  elevation  time  histories  respectively  for  fo  =  10.0 
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MHz.  The  most  prominent  difference  between  these  figures  and  the  earlier  ones  for  7.0  MHz 
(Figure  8a  and  b),  is  the  reduced  point-to-point  scattering  in  the  estimated  bearings  at  the  start  and 
end  of  the  run.  This  is  explained  by  the  narrowness  of  the  strip  wherein  scattering  centers  can 
contribute  to  the  signal.  The  intersection  of  this  strip  with  the  scattering  patch  at  the  start  and  end 
of  the  run  subtends  a  smaller  bearing  angle  at  the  receive  array,  giving  rise  to  the  smaller  spread  in 
estimated  bearings.  Another  difference  is  in  the  elevation  angles,  which  for  any  particular  time  are 
lower  than  those  found  with  a  7.0  MHz  critical  frequency.  These  elevation  angles  are  lower  since 
with  a  higher  critical  frequency  the  paths  do  not  penetrate  as  deeply  into  the  ionosphere. 


4.2.2  Effective  Scattering  Region 


From  the  previous  discussion  for  the  10.0  MHz  critical  frequency,  it  is  evident  that  the 
scattering  centers  which  are  capable  of  contributing  to  the  signal  are,  by  virtue  of  the  height  at 
which  they  intersect  the  propagation  path,  limited  to  a  narrow  region  which  is  centered  on  the 
perpendicular  bisector  of  the  great-circle  transmit-receive  path.  Equation  17,  which  gives  the 
propagation  distance  Z)  as  a  function  of  the  zenith  angle  <|)q,  was  used  to  calculate  the  boundaries 
of  this  region,  as  described  below. 

In  order  to  calculate  the  locations  of  these  boundaries,  a  computer  program  was  written 
which  takes  as  input  the  ionospheric  parameters:  height  of  maximum  h^,  half-width  y^,  and 
critical  frequency  /q;  the  transmitter  and  receiver  (X,  Y)  coordinates,  and  the  operating 
frequency  f. 


In  the  program,  the  (X,  Y)  coordinate  system  is  chosen  so  that  the  transmitter  and  receiver 
lie  along  the  y-axis,  with  the  origin  at  the  halfway  point  (See  Figure  12).  The  perpendicular 
bisector  to  the  great-circle  transmit-receive  path  is  then  the  X-axis,  i.e.,  the  equation  7=0.  The 
program  steps  the  X-value  over  the  range  of  values  desired  for  plotting.  The  corresponding  y- 
coordinate  of  the  lower  boundary  to  the  effective-scattering-region  strip  is  found  for  a  particular 
X,  by  a  search  procedure  which  involves  decreasing  Y  from  zero,  in  steps.  For  each  point  (X,  Y) , 
the  program  computes  the  geometric  distance  D'  and  compares  it  to  the  distance  D  found  using 
equation  (17).  The  horizontal  path  length  D'  expected  from  the  geometry  is  calculat^  from  the 
(X,  Y)  location  in  question,  together  with  the  transmit  and  receive  locations.  A  zenith  angle  is 
obtained  by  using  the  boundary  condition 


<j)Q  =  atan 


I  J 


(44) 


where  di  is  the  geometric  ground  distance  between  the  scatter  and  receive  locations.  The 
propagation  distance  D  corresponding  to  <1)q  is  then  computed  from  (17)  and  compared  with  the 
geometrically  obtained  value  D' .  Initially,  at  y  =  0,  D  will  be  greater  than  D  .  As  Y  becomes 
more  negative,  D  becomes  less,  and  eventually  drops  below  D  .  When  this  happens,  the  two  most 
recent  values  are  interpolated  to  determine  the  point  at  which  D  =  D' .  This  corresponds  to  the 
Tninitniim  value  The  maximum  value,  by  symmetry  is  given  by  ^ max  ^min’  SOme 

critical  frequencies  andX-values,  as  Y  is  made  more  negative,  the  ray  will  begin  to  penetrate  right 
through  the  ionosphere  before  D  (found  from  (44)  and  (17))  drops  as  low  as  the  geometric 
distance  D' .  This  situation  corresponds  to  the  geometric  distance  D' being  inside  the  skip  range  so 
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that  propagation  is  impossible.  The  penetration  condition  was  checked  for  at  each  step  in  the 
program,  by  invoking  equation  (19).  In  this  way,  the  boundaries  of  the  strip  of  allowed  scattering 
locations  are  found  for  plotting. 

Figure  12  is  the  resulting  plot  of  the  boundaries  of  the  strip  within  which  scattering  centers 
are  capable  of  contributing  to  the  signal  at  the  receiving  array.  The  minimum  and  maximum 
values  of  y  are  plotted  as  a  function  of  x,  for  several  values  of  critical  frequency  /q.  The  transmit- 
receive  geometry  and  ionospheric  parameters  are  the  same  as  used  previously. 


o 


'-750 


750 


X  COORDINATE  (km) 


Figure  12.  Boundaries  of  the  allowed  region  wherein  scattering  centers  can  contribute 
to  the  signal  seen  at  the  receiving  array,  for  various  critical  frequencies. 


Several  things  are  apparent  from  this  figure.  The  width  (i.e.,  y-extent)  of  the  allowed  region, 
when  the  critical  frequency  is  less  than  that  needed  to  support  great-circle  propagation  (e.g.,  6.0 
MHz),  increases  as  one  moves  closer  to  the  great  circle,  up  to  a  point  where  the  rays  begin  to 
penetrate  through  the  ionosphere  making  propagation  to  the  receiver  impossible.  There  is  also  a 
slight  widening  near  the  great  circle,  for  critical  frequencies  just  above  that  needed  to  support 
great-circle  propagation  (6.5  MHz).  However,  as  the  critical  frequency  increases  further,  an 
opposite  effect  occurs.  At  a  critical  frequency  of  10.0  MHz,  the  width  of  the  allowed  region 
decreases  slightly  as  the  great  circle  is  approached. 

The  width  of  the  allowed  region  where  it  crosses  the  great  circle  is  seen  from  Figure  12  to  be 
68  km,  for  a  critical  frequency  of  10.0  MHz.  This  is  in  good  agreement  with  the  67.5  km  width 
inferred  from  the  number  of  contributing  scattering  centers  (Section  4.2.1. 1). 


27 


4.23  Variation  of  Scattering  Efficiency  with  Height 

The  simulations  of  Section  4.2.1  which  were  performed  using  the  stronger  scattering  width 
height  dependence  given  by  equation  (29)  were  repeated  for  the  weaker  dependence  given  by 
equation  (30),  based  upon  an  excess  scatterer  electron  density  proportional  to  the  square  root  of 
the  ambient  electron  density  (5iV«  Jn).  The  results  were  found  to  be  essentially  the  same  as  for 
the  stronger  height  dependence  for  which  5N  oe  N.  However,  small  differences  existed. 

There  is  a  tendency  for  the  mean  elevation  of  the  direction  estimates  to  be  lower  when  the 
scattering  efficiency  varies  less  with  height  {8N  «>«  Jn),  which  is  expected  since  the  higher  angle 
rays  which  intersect  their  corresponding  scattering  center  higher  in  the  ionosphere  would  not  be 
as  strongly  favoured  with  this  height  dependence.  For  the  same  reason,  there  is  also  a  tendency  for 
the  rms  spread  in  the  elevation  angles  experienced  at  any  one  patch  position  to  be  slightly  greater 
for  the  reduced  height  dependence.  It  should  be  emphasized  that  these  differences  are  small,  so 
that  the  actual  form  of  the  scattering-width  height  dependence  plays  only  a  small  role  in 
determining  the  behaviour  of  the  field-aligned  scattering  model. 

Further  simulations  presented  in  this  paper  use  the  weaker  (equation  (30))  height 
dependence  for  the  effective  scattering  width. 


4.2.4  Multiple-Patch  Case 

In  any  real  high-latitude  propagation  situation  where  there  are  ionospheric  gradients  and 
small-scale  irregularities  present  as  a  result  of  the  gradient-drift  mechanism,  the  small-scale 
irregularities  will  be  grouped  into  specific  regions  or  patches  in  the  vicinity  of  the  gradients. 
Fiuther,  with  the  large-scale  irregularities  likely  to  exist  at  polar  latitudes  (e.g..  Figure  1)  there 
will  likely  be  more  than  one  region  of  high  gradient  present,  and  therefore  more  than  one 
scattering  patch  present.  The  number  of  patches  actually  observed  will  depend  on  their  position 
relative  to  the  transmitter  and  receiving  array,  as  well  as  the  propagation  conditions. 


As  an  initial  approach  to  multiple-patch  simulation  two  patches  are  considered,  each  200  km 
square,  with  similar  scattering  and  ionospheric  parameters.  The  ionospheric  parameters  assumed 
are  given  in  Section  4.1,  and  are  similar  to  those  used  for  the  single  patch  case  (Section  4.2.1). 

A  plan  view  of  the  geometry  used  for  many  of  the  simulations  is  given  in  Figure  13. 


Figure  13.  Plan  view 
of  modelled  two-patch 
geometry. 
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4.2.4.1  Initial  Simulations 


The  initial  simulations  used  a  critical  frequency  of  6.5  MHz,  just  above  the  minimum 
required  for  great-circle-path  propagation.  A  scattering-center  density  of  0.01  centers/km^  was 
used,  so  that  on  average  Aere  were  400  scattering  centers  in  each  patch. 

Figure  14  shows  the  time  history  of  the  estimated  plane-wave  direction,  described  in  terms 
of  bearing  and  elevation.  As  the  simulation  run  progressed,  the  observed  bearings  in  Figure  14  are 
directed  mostly  toward  the  nearer  patch,  first  directed  to  patch#  1,  then  to  patch  #2.  A  separate 
analysis  of  two  plane-wave  signals  shows  that  a  plane-wave  estimator  will  be  always  directed 
primarily  toward  the  stronger  signal,  provided  it  is  1  or  more  dB  above  the  weaker  one.  For  much 
of  the  time,  the  observed  bearing  jumps  between  the  two  patches,  suggesting  that  the  composite 
signals  from  the  patches  vary  randomly  in  amplitude  so  that  sometimes  one  patch  has  the  stronger 
signal,  sometimes  the  other.  The  corresponding  elevation  time  history  in  Figure  14  appears  to  be 
the  composite  of  two  arcs,  each  peaking  near  the  time  their  implied  patch  crosses  the  great  circle. 
Near  the  middle  of  the  run,  the  observed  elevations  appear  to  jump  between  the  implied  arc  for 
one  patch,  and  the  arc  for  the  other. 


Figure  14.  Time  histories  of  estimated  plane  wave  bearing  and  elevation, 
for  the  two-patch  geometry  of  Figure  13,  with  a  critical  frequency 
of  65  MHz  and  a  scattering-center  density  ofO.OlIkmr. 

It  may  be  thought  that  the  degree  of  randomness  in  the  composite  signal  coming  from  a 
single  patch,  and  therefore  the  uncertainty  as  to  which  patch  contributes  the  stronger  signal  thus 
capturing  the  plane-wave  direction,  is  a  result  of  the  small  number  of  scattering  centers  used  (an 
average  of  400  per  patch).  To  see  if  this  is  the  case,  the  scattering-center  density  was  increased 
tenfold,  to  0.1/km^,  (4000  scattering  centers  per  patch)  and  the  simulation  repeated.  No  dramatic 
change  was  seen  in  the  estimated  plane-wave  direction  histories  when  the  scattering-center 
density  was  increased  tenfold.  The  duration  of  time  where  both  patches  were  seen,  i.e.,  where  the 
estimated  plane-wave  direction  jumped  back  and  forth  between  the  general  directions  of  the  two 
patches  was  roughly  the  same  for  the  two  densities. 


29 


4.2.4.2  Dependence  on  Ionospheric  Electron  Density 

In  order  to  illustrate  the  strong  role  played  by  ionospheric  refraction  in  the  scattering  model, 
and  the  degree  of  complexity  that  can  result  in  a  direction-finding  record  when  more  than  one 
scattering  patch  is  present,  tests  were  run  using  various  critical  frequencies  both  above  and  below 
that  needed  to  support  great-circle  propagation.  The  geometry  of  Figure  13  was  used  for  these 
tests,  along  with  a  scattering-center  density  of  0.02/km^.  The  resultant  bearing  and  elevation 
histories  are  shown  in  Figure  15. 

For  a  critical  frequency  /q  of  5.0  MHz,  propagation  is  not  supported  for  scattering-center 
positions  within  440  km  of  the  great  circle,  represented  by  0°  bearing..  At  the  start  of  the  run, 
only  the  scattering  centers  in  the  furthest  patch  (#2)  contribute  to  the  signal,  and  when  this  patch 
moves  closer  than  440  km,  the  signal  disappears  altogether,  near  estimate  180.  Up  to  that  time,  the 
bearing  and  elevation  estimates  of  Figure  15  reflect  the  location  of  the  patch  #2  scatterers  which 
contribute  to  the  signal.  For  estimates  180  to  320,  there  is  no  propagation  for  scattering  centers 
from  either  patch,  and  thus  no  signal.  After  estimate  320,  patch  #1  moves  past  440  km  on  the 
other  side  of  the  great  circle,  and  scattering  centers  from  that  path  begin  to  contribute  to  a  signal. 
The  plane-wave  estimates  for  the  signal  reflect  the  scattering  positions  within  patch  #1  that  are 
more  than  440  km  from  the  great  circle.  Thus,  for  a  5.0  MHz  critical  frequency,  only  the  further  of 
the  two  patches  is  seen,  and  for  only  part  of  the  time,  owing  to  the  large  skip  distance  imposed  by 
the  ionosphere. 

At  /o  =  6.0  MHz,  propagation  is  not  supported  for  scattering  positions  within  100  km  of  the 
great  circle.  At  the  beginning  of  the  6.0  MHz  run,  signals  from  both  patches  are  present,  and  the 
resultant  bearings  and  elevations  of  Figure  15  are  directed  toward  the  nearer  patch  (#1).  After  the 
first  80  estimates,  the  patches  move  so  that  some  scattering  centers  inside  patch  #1  are  closer  than 
100  km  from  the  great  circle  and  do  not  contribute  to  its  signal.  At  this  point  the  estimated  plane- 
wave  signal  direction  switches  to  patch#2  which  now  has  the  stronger  contribution.  The  estimated 
directions  continue  to  follow  the  locations  of  those  scattering  centers  inside  patch  2  which  remain 
outside  100  km  of  the  great  circle.  Patch  #1  crosses  the  great  circle  and  begins  to  move  outside  the 
100-km  limit  on  the  other  side  of  the  great  circle.  Patch  #1  begins  to  again  contribute  to  the  signal, 
while  patch  #2  which  is  now  moving  inside  the  100-km  limit  on  the  other  side  of  the  great-circle 
begins  to  have  its  contribution  reduced.  For  a  period  (estimates  220  to  280),  the  contributions  of 
both  patches  are  fairly  close  in  mean  power,  so  that  random  fluctuations  cause  the  resultant  plane- 
wave  direction  to  jump  back  and  forth  between  the  positions  of  the  contributing  scattering  centers 
in  each  patch.  After  estimate  280,  patch  #2  moves  completely  inside  the  100-km  limit  and  patch 
#1,  outside  the  100-km  limit  on  the  other  side.  The  direction  estimates  reflect  the  location  of  patch 
#1  at  this  point.  This  continues  until  patch  #2  crosses  the  great  circle  and  moves  outside  the  100- 
km  limit,  whereupon  it,  as  the  nearer  of  the  two  patches,  begins  to  control  the  estimated  plane- 
wave  direction  (after  estimate  410). 


At  /o  =  6.5  MHz,  propagation  is  supported  for  all  scattering  locations.  The  estimated  plane- 
wave  directions  are  controlled  by  the  patch  that  has  the  stronger  composite  signal.  This  depends 
not  only  on  the  number  of  scattering  centers  in  the  patch,  but  also  the  effective  scattering  width 
for  each  scattering  center  (and  hence  the  height  of  scattering),  and  the  ionospheric  focusing  that 
occurs  for  that  path.  In  addition,  the  geometric  path  loss  involved  also  favours  the  shorter  paths. 
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Figure  15.  Estimated  plane-wave  bearing  and  elevation  time  histories  for  the  two-patch 
geometry  of  Figure  13,  for  various  critical  frequencies,  and  an  8.0  MHz  transmitted  frequency. 


At  the  start  and  end  of  the  run  illustrated  in  Figure  15,  the  estimated  plane- wave  directions  are 
directed  toward  the  closer  of  the  two  patches.  There  is  a  long  period  in  the  middle  of  the  run 
where  the  estimate  direction  moves  between  the  two  patches,  due  to  random  fluctuations  in  the 
roughly  equal  signal  powers  now  favouring  one  patch,  then  the  other. 


As  /o  is  increased  to  7.0  and  then  10.0  MHz  in  Figure  15,  the  observed  directions  continue 
to  favour  the  closer  patch.  However,  the  range  of  times  over  which  both  patches  can  be  seen 
increases,  suggesting  that  the  composite  powers  fiom  each  patch  are  more  nearly  equal  for  a 
wider  range  of  positions  than  previously.  The  main  difference  between  the  results  for  these  critical 
frequencies  and  those  of  the  6.5  MHz  critical  frequency  is  the  absence  of  those  the  strong 
ionospheric  focusing  that  occurs  for  the  midpoint  paths  near  great  circle  whose  ground  distance  is 
close  to  that  of  the  skip  distance.  This  weakening  of  the  ionospheric  focusing  effect  clearly 
reduces  the  dependence  of  the  composite  signal  strength  upon  patch  location. 


4.2  J  Effect  of  Different  Geometries 


The  geometries  modelled  so  far  have  been  restricted  to  patches  moving  along  the 
perpendicular  bisector  of  the  great-circle  path.  This  is  generally  not  the  case,  and,  as  shown  in 
Section  4.2.2,  the  perpendicular  bisector  occupies  a  privileged  position  at  the  midpoint  of  the 
allowed  scattering  region.  It  is  therefore  enlightening  to  consider  other  geometries. 


Figure  16  shows  a  modelled  geometry  where  two  scattering  patches  follow  one  another 
across  the  midpoint  of  the  great-circle  path,  in  a  direction  that  is  26.5°  off  the  perpendicular 
bisector,  corresponding  to  a  2-km  shift  in  x  and  1-km  shift  in  y  between  estimates.  A  critical 
frequency  of  10.0  MHz  was  used  for  this  test,  along  with  a  scattering-center  density  of  0.02/km  . 
The  resultant  bearing  and  elevation  history  is  shown  in  Figure  17. 


Figure  16.  Plan  view  of 
two-patch  geometry  used 
in  obtaining  Figure  17. 


km 


When  the  bearing  time  history  of  Figure  17  is  compared  with  that  of  Figure  15  for  a  critical 
ftequency  /o  of  10.0  MHz,  large  differences  are  seen  although  the  change  in  geometry  is  not 
great.  The  bearing  estimates  in  Figure  17  switch  quite  abruptly  from  one  patch  to  the  other  near 
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Figure  17.  Time  history 
of  estimated  plane-wave 
bearing  for  fo  =  10  MHz, 
for  the  geometry  of 
Figure  16. 
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the  midpoint  of  the  run,  while  those  in  Figure  15  jump  between  patches,  allowing  both  patches  to 
be  seen  at  the  same  time  over  most  of  the  run.  The  abrupt  change  in  Figure  17  can  be  related  to  the 
scattering  requirement  that  the  scattering-center  locations  in  order  to  intersect  the  propagation 
paths  inside  Ae  ionosphere  must  lie  within  a  narrow  strip  centered  on  the  perpendicular  bisector 
of  the  great-circle  path,  as  illustrated  in  Figure  12.  Thus,  for  the  geometry  of  Figure  16,  each  patch 
is  seen  in  turn  as  it  intersects  the  strip,  and  only  for  a  brief  time  near  the  middle  of  the  run  are  both 
patches  seen  simultaneously. 

As  further  influence  of  the  effect  of  a  changing  drift-motion  direction,  the  situation  of  two 
patches  following  one  another  along  a  patch  crossing  the  midpoint  of  the  great-circle  path  at 
63.5®  off  the  perpendicular  was  modelled,  corresponding  to  a  1  km  shift  in  x  and  2  km  shift  in  y 
between  estimates.  A  lower  critical  frequency  /q,  of  7.0  MHz  was  used.  For  this  frequency,  Figure 
12  shows  a  wider  strip  of  permissible  scattering  locations  as  the  rays  penetrate  further  into  the 
ionosphere.  Figure  18  illustrates  the  geometry,  and  Figure  19  shows  the  resultant  time  history  of 
estimated  plane- wave  bearings. 


Figure  18.  Plan  view 
of  two-patch  geometry 
used  in  Figure  19. 
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As  in  the  previous  case  the  estimated  directions  follow  one  patch  then  the  other,  as  they  pass 
over  the  strip  of  allowed  scattering  locations.  As  seen  in  Figure  19,  there  is  a  time  in  the  middle  of 
the  run  when  the  two  patches  are  on  either  side  of  the  strip  so  that  neither  patch  is  seen,  and  at  the 
start  and  end  of  the  run  the  two  patches  cannot  be  seen,  as  they  are  to  one  side  of  the  strip. 

4.3  Summary  and  Discussion  of  Findings 

The  scattering  model  was  found  to  produce  a  composite  received  signal  at  the  array  that  fit  a 
plane  wave  surprisingly  well  most  of  the  time,  although  the  constituent  signal  contributions 
arrived  from  many  directions  from  different  scattering  centers.  When  the  model  is  run  a  number 
of  times  under  the  same  conditions,  the  estimated  plane-wave  directions  are  distributed  in 
direction  so  as  to  reflect  the  spread  in  direction  of  the  constituent  signals  from  the  contributing 
scattering  centers. 

In  general,  the  spread  in  bearing  for  a  square  patch  of  scattering  centers  is  much  larger  than 
the  spread  in  elevation.  The  spreads  depend  not  only  on  the  actual  positions  of  the  scatterers  in  the 
patch  relative  to  the  receive  array,  but  also  on  the  propagation  conditions  which  determine  the 
height  of  scattering,  and  the  scattering  locations  for  which  propagation  is  possible  (Figures  10, 
11).  The  observations  are  not  strongly  dependent  on  the  form  of  the  height  dependence  of  the 
scattering  cross  section. 

For  the  horizontally  stratified  ionosphere  modelled,  the  scattering  centers  that  can  contribute 
to  the  received  signal  are  restricted  to  a  strip  of  varying  width,  centered  about  the  perpendicular 
bisector  of  the  great-circle  transmit-receive  path  (Figure  12).  The  width  of  this  snip  for 
ionospheric  densities  well  above  that  needed  to  support  great-circle  propagation  is  quite  narrow, 
less  than  the  dimensions  anticipated  for  many  patches  of  scatterers.  As  the  ionospheric  electron 
density  becomes  less,  the  strip  widens,  until  the  point  is  reached  where  great-circle  propagation  is 
no  longer  supported,  whereupon  it  splits  in  two.  The  scatterers  contributing  to  the  signal  at  any 
time  are  located  in  the  intersection  of  the  patches  present  with  the  allowed  regions  (as  defined  by 
the  ionospheric  conditions)  for  the  scatterers  to  contribute.  For  a  non-horizontally  stratified 
ionosphere,  the  shape  of  the  allowed  region  will  be  more  complex  than  a  simple  strip,  and  will 
depend  strongly  on  the  ionospheric  contouring  present. 
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In  the  case  of  two  or  more  patches  of  scatterers,  the  patch  which  has  the  strongest  composite 
signal  at  the  receiving  array  will  control  the  estimated  plane-wave  direction.  The  present  model 
assumes  that  the  scatterers  are  regenerated  randomly  between  consecutive  estimates,  based  upon 
the  short  lifetimes  expected  for  small  field-aligned  irregularities.  As  a  result,  the  resultant 
composite  signal  for  any  one  patch  will  fluctuate  substantially  in  its  strength,  from  estimate  to 
estimate.  When  conditions  are  such  that  two  patches  produce  signals  which  are  fairly  close  in 
their  mean  power,  the  estimated  plane-wave  direction  will  jump  from  one  patch  to  the  other. 
Depending  on  the  geometry  and  propagation  conditions,  two  or  more  patches  may  be  seen 
together  in  this  way  for  some  time  (Figure  15, 7  and  10  MHz). 

When  the  ionospheric  electron  density  is  lower  than  that  needed  for  great-circle 
propagation,  scatterers  along  the  great-circle  will  not  contribute  to  the  observed  signal,  but  others 
further  away  will  contribute,  provided  they  lie  inside  the  allowed  region.  Patches  of  scatterers 
moving  through  the  transmitter  receiver  path  will  thus  cause  the  time  history  of  received  signal 
directions  to  be  quite  complex,  and  highly  dependent  of  the  ionospheric  propagation  cutoff  and 
the  geometry  of  the  transmitter  and  receiver  with  respect  to  the  patches  and  their  motion  (Figure 
15, 5  and  6  MHz,  Figures  17  and  19). 
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5.0  TESTS  OF  THE  MUSIC  ALGORITHM  APPLIED  TO  THE 

FIELD-ALIGNED  SCATTER  MODEL 

In  considering  the  implications  of  aspect-sensitive  scattering  upon  direction-finding,  the 
propagation  model  and  the  DF  technique  cannot  be  viewed  in  isolation,  but  need  to  be  treated 
together.  In  the  previous  section,  the  implications  of  the  aspect-sensitive  scattering  model  for 
direction-finding  were  reviewed,  for  the  case  of  a  plane-wave  fitting  (i.e.,  interferometric)  DF 
technique.  Those  results  are  useful  in  establishing  not  only  the  effects  of  the  propagation  model 
(i.e.  scattering  and  ionospheric  refraction)  but  also  the  way  these  effects  are  seen  by  the  DF 
system,  when  a  plane-wave  technique  is  employed.  Such  results  lend  themselves  to  comparisons 
with  real  data  obtained  by  systems  using  interferometric  techniques,  such  as  that  gathered  by  the 
Kestrel  experiment  [9,10]. 

In  the  present  section,  the  implications  of  scattering  upon  direction  finding  are  extended  to  a 
class  of  direction-finding  algorithms  of  which  the  MUSIC  algorithm  is  representative  These 
algorithms  are  coming  into  use  in  sampled-aperture  systems,  and  are  capable  of  obtaining  a 
number  of  signal  directions  simultaneously.  The  signal  directions  can  represent  different 
transmissions,  or  a  number  of  coexisting  propagation  paths  for  the  same  transmission.  These 
algorithms  consider  a  small  number  of  directions,  less  than  the  number  of  antennas  in  the 
receiving  array,  so  that  the  number  of  propagation  paths  which  can  be  exactly  identified  is  limited. 
In  the  case  of  scattering,  we  arc  dealing  with  a  large  number  of  paths,  many  more  than  can  be 
specified.  How  the  MUSIC  algorithm  will  respond  to  this  situation  needs  to  be  determined. 


5.1  Single  Scattering  Patch 

5.1.1  Test  Parameters 

Initial  tests  with  the  MUSIC  algorithm  deal  with  a  single  scattering  patch.  The  geometry  and 
ionospheric  parameters  are  the  same  as  those  previously  used  (Sections  4.1, 4.2.1,  and  Figure  7). 
The  effective  scattering-width  dependence  on  height  is  assumed  to  be  proportional  to  the 
ionospheric  electron  density,  as  given  by  equation  (30). 

The  MUSIC  algorithm  requires  a  series  of  signal  samples  to  be  taken  for  each  direction 
estimate.  In  the  present  tests,  50  samples  were  generated  for  this  purpose,  for  each  estimate.  The 
samples  were  taken  2  ms  apart,  and  a  20-Hz  Doppler  spread  was  assumed  for  the  signal  from  each 
scattering  center  (Section  2.6.3). 

Consecutive  direction  estimates  were  made  with  a  4  km  shift  in  the  patch  position  between 
estimates.  This  is  equivalent  to  an  estimate  every  8  seconds  for  the  500m/sec  patch  motion  typical 
of  polar-cap  ionospheric  convection. 
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5.1.2  Simulation  Results 


Figure  20  presents  the  simulation  results.  Time  histories  of  the  estimated  bearings  and 
elevations,  found  for  the  single-patch  case  using  the  MUSIC  algorithm,  are  presented,  for  several 
settings  of  the  MUSIC  algorithm,  i.e.,  different  numbers  of  assumed  signal  directions.  These 
results  may  be  compared  with  those  of  Figure  8a  and  b  of  Section  4.2.1,  obtained  for  the  same 
geometry  and  ionospheric  conditions  using  plane-wave  fitting. 

The  results  found  with  MUSIC  depend  strongly  on  the  number  of  assumed  signal  directions. 

From  Figure  20,  when  a  single  direction  is  assumed,  the  resultant  direction  estimates  tend  to 
follow  the  central  portion  of  the  patch.  There  is  some  random  scatter  in  the  point-to-point 
estimates,  but  not  so  much  that  the  directions  can  be  said  to  come  from  anywhere  in  the  patch. 
The  spread  in  directions  is  much  less  than  the  angle  subtended  by  the  patch,  in  contrast  to  the 
result  found  with  the  plane-wave  fitting  procedure  in  Figures  8a  and  b,  where  the  spread  was  more 
representative  of  patch  size. 

When  two  signal  directions  are  assumed,  the  estimated  directions  of  Figure  20  tend  to 
represent  the  two  sides  of  the  patch  as  seen  at  the  receiving  array.  There  are,  in  effect  two  traces, 
on  the  bearing  and  elevation  time  history  plots,  representing  the  two  sides.  At  the  start  of  the  run, 
the  lower  trace  on  the  bearing  plot  represents  the  side  of  the  patch  nearest  to  the  great-circle,  and 
corresponds  to  the  elevation  trace  with  the  higher  elevation.  Each  side  is  associated  with  a  bearing 
trace  that  becomes  lower  with  time,  and  an  elevation  trace  that  increases  up  to  the  point  where  the 
side  crosses  the  great-circle  transmit-receive  path  after  which  it  decreases.  The  separation  of  the 
two  bearing  traces  is  representative  of  the  angle  subtended  by  the  patch  at  the  receiver.  However 
the  two  traces  would  appear  to  suggest  two  smaller  patches  rather  than  a  single  extended  patch. 

In  addition  to  the  two  traces  representing  the  edges  of  the  patch,  there  are  several  shorter 
traces  that  appear  during  the  run,  which  bear  no  clear  relationship  to  the  patch  position,  and  which 
have  relatively  fixed  azimuths.  There  are  two  likely  causes  for  these  extra  traces.  One  possibility 
is  due  to  the  large  sidelobes  which  are  beginning  to  appear  in  the  receiving  array  pattern  at  the  8 
MHz  operating  frequency.  Directions  of  sidelobes  whose  amplitude  whose  amplitude  is  close  to 
that  of  the  main  lobe  may  be  identified  as  signal  directions  whenever  the  actual  si^d  comng 
from  a  large  number  of  other  directions  is  sufficiently  nonplanar.  The  second  possibility  is  an 
artifact  of  the  search  procedure  employed  for  finding  the  signal  directions,  which  may 
occasionally  cause  certain  boundary  points  of  the  direction  cosines  searched  to  be  identified  as 
signal  directions. 

When  the  number  of  assumed  signal  directions  is  increased  to  three  and  then  to  four,  Ae 
direction  estimates  tend  to  fill  in  the  area  between  the  two  traces  found  for  the  two-direction 
estimation.  The  fact  that  the  signals  come  fi-om  a  single  extended  patch  becomes  fairly  clear.  The 
direction  estimates  cover  a  region  whose  extent  is  representative  of  the  angle  subtended  by  the 
patch  at  the  receiving  array.  Further,  the  subtended  angle  is  much  better  defined  in  these  cases 
than  it  is  for  the  plane-wave-fitting  result  in  Figures  8a  and  b.  At  the  same  time,  the  extra  traces 
possibly  arising  from  array  sidelobe  effects  and  algorithm  artifacts  become  much  stronger, 
enough  so  that  measures  to  remove  them  should  be  considered. 
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1  signal  direction  assumed 


3  signal  directions  assumed 


4  signal  directions  assumed 


Figure  20.  Time  histories  of  estimated  bearing  and  elevation,  for  the  single-patch 
geometry  of  Figure  7,  and  a  critical  frequency  of  7.0  MHz,  as  determined  by  the 
MUSIC  algorithm,  for  various  numbers  of  assumed  signal  directions. 
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5.2  Two  Scattering  Patches 


5.2.1  Test  Parameters 

The  geometry  used  for  the  MUSIC  simulations  is  the  same  as  that  used  with  the  plane-wave 
fitting,  shown  in  Figure  13.  The  transmitter  location  and  receiving  array  parameters  are  the  same 
as  for  the  single-patch  simulations,  i.e.,  similar  to  those  of  the  Kestrel  measurements.  The 
frequency  of  operation  is  8.0  MHz. 

Apart  from  the  two-patch  geometry,  all  other  parameters  are  the  same  as  those  used  for  the 
single-patch  case  (Section  5.1.1).  The  critical  frequency  is  7.0  MHz,  which  is  slightly  higher  than 
the  minimum  needed  to  support  propagation  over  the  great-circle  transmit-receive  path. 


5.2.2  Test  Results 

The  simulations  were  run  with  for  the  assumed  number  of  signal  directions  set  to  2, 3, 4,  and 
5.  Figure  21  gives  the  time  histories  of  the  bearing  and  elevation  estimates  obtained.  These  results 
may  be  compared  with  those  found  using  plane- wave  fitting  with  the  same  geometry  and  (7.0 
MHz)  critical  frequency,  in  Section  4.2.4. 1,  Figure  14. 

The  main  difference  between  the  present  results  and  those  of  plane-wave  fitting,  is  that  both 
patches  are  seen  at  all  times  during  the  run  using  the  MUSIC  algorithm  with  more  one 
assumed  direction.  The  plane-wave-fitting  technique  tends  to  favour  the  nearer  patch,  giving  the 
farther  patch  limited  visibility  especially  at  the  start  and  end  of  the  run. 

When  two  signal  directions  are  assumed  by  the  MUSIC  algorithm.  Figure  21  shows  the 
resultant  signal  directions  to  be  slightly  scattered,  about  two  traces  centered  on  the  directions  to 
the  centers  of  the  two  patches.  The  spreading,  as  in  the  single-patch  case  for  one  assumed 
direction,  is  less  than  that  observed  with  the  plane-wave  fitting. 

When  the  number  of  assumed  signal  directions  increases  above  the  number  of  patches,  the 
estimated  directions  become  more  spread,  so  that  they  cover  the  angles  subtended  by  the  two 
patches.  For  four  assumed  directions,  the  estimates  obtained  tend  to  follow  the  edges  of  the  two 
patches  as  seen  at  the  receiving  array,  although  this  effect  is  not  as  pronounced  as  for  a  single 
patch  and  two  assumed  directions.  For  five  assumed  directions,  the  estimates  tend  to  fill  in  the 
angles  subtended  by  the  patches,  so  that  reasonable  estimates  of  patch  size  may  be  obtained. 

As  in  the  single-patch  case,  when  the  number  of  assumed  directions  is  increased  beyond  the 
number  of  patches,  additional  traces  begin  to  appear  which  bear  little  relation  to  patch  position. 
The  bearings  of  these  traces  appear  to  change  less  with  time  than  those  associated  with  patch 
position. 
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2  Signal  directions  assumed 


3  signal  directions  assumed 


4  signal  directions  assumed 


5  signal  directions  assumed 


Figure  21.  Time  histories  of  estimated  bearing  and  elevation,  for  the  two-patch 
geometry  of  Figure  13,  and  a  critical  frequency  of  7.0  MHz,  as  determined  by 
the  MUSIC  algorithm,  for  various  numbers  of  assumed  signal  directions. 
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5.3  Summary  and  Discussion  of  Findings 


The  simulations  show  that  the  MUSIC  algorithm,  when  applied  to  scatter  propagation  from 
an  extended  source,  has  certain  advantages  over  the  plane-wave  fitting  algorithm.  When  several 
patches  are  present  at  the  same  time,  it  is  easier  to  detect  the  presence  of  each  with  the  MUSIC 
algorithm  which  can  look  for  several  signal  directions  at  once.  Further,  a  clearer  picture  of  the 
presence  and  extent  of  an  extended  scattering  region  is  obtained  with  MUSIC  in  a  relatively  short 
sampling  time. 

However,  the  MUSIC  algorithm  appears  to  be  prone  to  providing  false  direction  estimates, 
exemplified  by  the  additional  traces  of  Figures  20  and  21.  These  additional  traces  are  most  likely 
a  result  of  the  mismatch  between  the  small  number  of  signal  directions  assumed  by  MUSIC  and 
the  very  large  number  of  actual  signal  directions  from  the  extended  scattering  patches.  As  will  be 
shown  in  the  next  section,  they  are  more  likely  to  appear  when  the  receiving  array  sidelobes  are 
large.  This  implies  that  they  are  due  to  a  sidelobe  direction  sometimes  being  identified  as  a  true 
signal  direction,  which  is  more  likely  to  occur  when  the  sidelobe  gain  is  close  to  that  of  the  main 
lobe. 
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6.0  TESTS  OF  THE  MUSIC  AGORITHM  APPLIED  TO  FIELD-ALIGNED 

SCATTERING  AND  SPORADIC-E 

In  the  Kestrel  tests  conducted  in  the  polar  cap  region  [9],  sporadic  E  was  identified  as  the 
means  of  propagation  over  the  Thule- Alert  path  slightly  more  than  1  percent  of  the  time.  This 
propagation  mode  was  characterized  by  very  small  bearing  departures  from  great-circle,  in 
contrast  to  the  other  propagation  modes  which  were  largely  associated  with  scattering.  Therefore 
it  is  highly  desirable  that  this  mode  be  identified  when  present,  and  used  to  obtain  emitter 
bearings. 

As  mentioned  in  Section  4.2.4. 1,  a  separate  analysis  of  the  performance  of  plane- wave 
fitting  in  the  presence  of  two  signals  indicates  that  the  plane-wave  fitting  technique  always  yields 
a  direction  close  to  that  of  the  stronger  signal,  provided  the  stronger  signal  is  1  or  more  dB  above 
the  weaker  signal.  Thus,  when  several  propagation  modes  were  present,  the  Kestrel  equipment 
which  used  plane-wave  fitting  always  selected  the  stronger  signal.  Sudden  jumps  in  the  signal 
direction  were  a  frequent  occurrence  in  the  Kestrel  data;  these  were  likely  due  to  changes  in  the 
relative  strength  of  different  propagation  modes.  The  nature  of  many  of  these  signals  suggests  that 
they  were  likely  due  to  scattering  from  extended  patches. 

The  plane-wave  fitting  restricted  visibility  of  sporadic-E  to  those  times  when  it  was  the 
strongest  mode.  It  is  worth  examining  other  techniques  such  as  MUSIC  that  may  increase  the 
visibility  of  this  useful  mode  and  therefore  the  percentage  of  time  that  it  may  be  used.  The 
simulations  presented  in  this  section  examine  sporadic-E  propagation  along  with  field-aligned 
scattering  propagation  of  greater  strength,  as  seen  by  a  receiving  array  employing  the  MUSIC 
algorithm. 


6.1  Test  Parameters 

The  path  geometry  was  similar  to  that  used  for  the  scatter  simulations  of  the  previous 
section,  which  was  based  on  the  Kestrel  measurements.  The  initial  receiving  array  geometry  was 
likewise  similar  to  that  of  the  Kestrel  array,  i.e.,  7  antennas,  selected  from  24  antennas  uniformly 
placed  around  a  circle  of  25  m  radius;  However,  later  simulations  used  other  geometries.  The 
operating  frequency  was  8  MHz. 

The  ionospheric,  scatterer,  and  signal  sampling  parameters  were  the  same  as  in  Section  5. 
The  spoiadic-E  reflection  coefficient  was  initially  set  so  that  the  received  signal  power  from 
sporadic  E  was  7  to  10  dB  below  the  power  received  from  field-aligned  scattering. 

In  general,  the  simulations  were  repeated  with  different  numbers  of  assumed  signal 
directions  to  be  estimated. 
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6.2  Test  Results 

6.2.1  Single  Scattering  Patch,  Sporadic  E  and  Kestrel  Array  Geometry 


Figure  22.  Plan  view  of  geometry  used  for  single-patch-plus-sporadic  E  simulations. 

A  plan  view  of  the  geometry  used  these  initial  simulations  is  shown  in  Figure  22.  The 
ionospheric  critical  frequency  was  7.0  MHz,  which  allowed  F-layer  propagation  along  the  great 
circle.  A  sporadic-E  reflection  coefficient  of  0.001  was  used;  this  value  provided  a  sporadic-E 
signal  power  that  was  7  dB  below  the  mean  signal  power  from  the  scattering  patch.  A  series  of 
250  estimates  were  made,  with  a  2-km  shift  in  patch  position  between  estimates. 

The  results  of  the  MUSIC  analysis  are  shown  in  Figure  23.  Hme  histories  of  the  estimated 
bearings  and  elevations  are  shown,  for  two,  three,  and  four  assumed  signal  directions. 

When  only  two  signal  directions  are  assumed,  one  of  those  directions  follows  the  moving 
scattering  patch,  and  the  other  fluctuates  about  the  (0°  bearing,  16°  elevation)  direction  of  the 
sporadic-E  signal.  There  is  a  large  amount  of  scatter  in  both  traces,  so  that  even  the  sporadic-E 
trace  does  not  provide  an  accurate  bearing  over  the  short  term.  (A  few  of  the  low-elevation  points 
appear  to  lie  exactly  along  a  12°  elevation  line;  these  points  correspond  to  bearings  near  0°  and 
may  be  an  artifact  of  the  direction-assignment  procedure.) 

When  the  number  of  signal  directions  in  each  estimate  is  increased  to  three,  the  directions 
begin  to  be  better  defined.  Two  of  the  traces  appear  to  follow  the  edges  of  the  patch,  as  previously 
noted  in  Figure  20  Section  5.1.2  for  the  single  scattering  patch  only  with  two  estimated  directions. 
The  other  trace  remains  fixed  on  the  sporadic-E  direction,  providing  a  good  estimate  of  the 
transmitter  bearing.  However  near  the  start  and  end  of  the  run,  when  the  patch  is  well  away  from 
the  sporadic-E  bearing,  there  is  a  tendency  for  the  third  direction  estimate  to  follow  the  scattering 
patch  rather  than  the  weaker  sporadic-E  signal.  The  sporadic-E  trace  only  becomes  clearly  visible 
when  the  scattering  patch  moves  closer  in  direction. 
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2  Signal  directions  assumed 


4  signal  directions  assumed 


Figure  23.  Time  histories  of  estimated  bearing  and  elevation,  for  the  geometry 
of  Figure  22,  and  a  critical  frequency  of  7.0  MHz,  as  determined  by  the  MUSIC 
algorithm,  for  various  numbers  of  assumed  signal  directions. 


When  four  signal  directions  are  assumed,  the  sporadic-E  trace  is  visible  throughout  the  run 
and  very  constant  in  azimuth.  The  scattering  patch  uses  the  remaining  three  directions  estimated, 
so  that  it  appears  in  the  time  histories  as  a  range  of  bearing  and  elevations  which  reflect  the  actual 
extent  of  the  patch  subtended  at  the  receiving  array.  Near  the  center  of  the  run,  an  additional  trace 
appears,  which  is  likely  an  array  sidelobe  effect  resulting  from  signals  amving  near  0°  bearing 
falling  in  a  sidelobe  of  the  array  steered  to  45°  bearing,  and  being  incorrectly  attributed  to  that 
direction.  (The  calculated  radiation  patterns  for  the  Kestrel  array  geometry  at  8  MHz  show  that 
when  the  array  is  steered  to  45  °  bearing  and  a  low  elevation  angle,  a  large  sidelobe  exists  near  0 
bearing  and  45°  elevation,  which  is  the  scattering  patch  direction  near  the  center  of  the  run). 
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6.2.2  Two  Scattering  Patches,  Sporadic  E  and  Kestrel  Array  Geometry 


Figure  24.  Plan  view  of  geometry  used  for  two-patch-plus-sporadic  E  simulations. 


Figure  24  shows  the  geometry  used  for  many  of  the  two-patch  simulations.  In  this  geometry, 
the  patches  follow  one  another  across  the  midpoint  of  the  great-circle  transmit-receive  path,  at 
right  angles  to  the  great  circle.  The  sporadic  E  reflection  coefficient  of  0.001  resulted  in  the 
sporadic-E  signal  being  approximately  10  dB  below  the  signal  received  from  the  two  scattering 
patches.  A  simulation  run  consisted  of  250  estimates,  with  a  4-km  shift  in  the  patch  location 
between  estimates. 

The  resulting  time  histories  of  estimated  bearing  and  elevation  are  shown  in  Figure  25  for 
two,  three,  four,  and  five  assumed  signal  directions. 

When  two  signal  directions  are  assumed  by  the  MUSIC  algorithm,  the  estimates  tend  to 
follow  the  positions  of  the  two  scattering  patches.  For  the  main  part  of  the  run,  the  sporadic-E 
propagation  is  masked  by  the  stronger  signals  from  the  scattering  patches,  with  one  assumed 
direction  allocated  to  each  of  the  scattering  patches.  Near  the  beginning  and  end  of  the  run, 
however,  the  signal  from  the  furthest  scattering  patch  is  weaker,  and  the  sporadic-E  propagation 
(near  0°  bearing)  begins  to  be  visible. 

When  three  signal  directions  are  assumed,  the  additional  direction  is  allocated  to  the 
sporadic-E  propagation  so  that  it  becomes  visible  throughout  the  run.  However  the  directions 
observed  are  somewhat  scattered  about  the  true  direction  of  0°  bearing.  When  the  assumed 
directions  are  increased  to  four,  the  sporadic-E  propagation  becomes  more  clearly  defined,  but 
some  scattering  in  diitction  remains;  two  of  the  assumed  directions  tend  to  be  appropriated  by  the 
nearer  scattering  patch  and  one  by  the  further  patch.  With  five  assumed  directions  both  patches  are 
well-defined  by  two  directions;  the  remaining  direction  is  then  free  to  provide  a  fairly  well- 
defined  sporadic-E  trace  with  little  scattering  in  bearing. 


45 


BEARING  (dagrMs)  _  BEARING  (dagrtas)  BEARING  j;dagraas)  BEARING  j;dagraas) 


2  Signal  directions  assumed 


4  signal  directions  assumed 


»1 

m  ' 
• 

•  , 

b. 

■gs- 

y 

UJ  * 

T. • .  ■ . ••  ' . .* 

II  '  ab  lAo  i4o  200  2. 

» 

^  '  ab  '  ido  2do  lio 

5  signal  directions  assumed 


Figure  25.  Time  histories  of  estimated  bearing  and  elevation,  for  the  two-patch  plus 
geometry  of  Figure  24,  and  a  critical  frequency  of  7.0  MHz,  as  determined  by  the  MUSIC 
algorithm  using  the  Kestrel  array,  for  various  numbers  of  assumed  signal  directions. 
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6.23  Effect  of  Antenna  Array  Geometry 

In  order  to  explore  the  impact  of  geometry  on  DF  performance  in  the  presence  of  scattering 
and  sporadic-E  propagation,  other  array  geometries  were  modelled.  The  log-spiral  geometry  of 
the  Vortex  array  currently  being  used  for  sampled-aperture  DF  tests  in  the  Ottawa  area  was  used. 

The  Vortex  array  consists  of  fifteen  vertical  whip  elements,  with  the  first  element  at  a  radius 
of  14.76m  from  the  spiral  center.  The  position  of  each  succeeding  element  is  increased  in  radius 
by  1.23  times  and  rotated  clockwise  about  the  center  by  an  angle  of  27c/7.  Normally  a  consecutive 
sequence  of  elements  (less  than  fifteen)  is  chosen  for  a  particular  period  of  operation,  with  the 
sequence  position  within  the  array  determined  by  the  operating  frequency,  inner  elements  being 
selected  for  high  frequencies  and  outer  elements  for  low  frequencies. 


5  signal  directions  assumed 


6  signal  directions  assumed 


7  signal  directions  assumed 


Figure  26.  Time  histories  of  estimated  bearing  and  elevation,  for  the  geometry  of  Figure  24 
and  a  critical  frequency  of  7.0  MHz,  as  determined  by  the  MUSIC  algorithm  using  the 
10  middle  Vortex  array  elements,  for  various  numbers  of  assumed  signal  directions. 
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As  an  initial  selection,  ten  consecutive  elements  from  the  middle  portion  of  the  fifteen- 
element  array  were  selected  for  study.  This  corresponded  to  an  array  for  which  the  smallest  radius 
was  27.47  m  and  the  largest,  176.99  m.  These  dimensions  and  most  of  the  resultant  element 
spacings  are  considerably  larger  than  those  of  the  Kestrel  array. 

Figure  26  shows  the  time  histories  for  the  bearing  and  elevation  estimated  by  MUSIC  with 
5,  6,  and  7  assumed  signal  directions,  for  the  10-element  Vortex  array  geometry.  The  propagation 
geometry  and  parameters  used  were  those  of  the  two-patch  case  covered  in  the  previous  section. 

The  capacity  of  this  array  for  seeing  the  weaker  sporadic-E  propagation  path  in  the  presence 
of  the  relatively  stronger  scattering  propagation  is  disappointing.  The  sporadic-E  path  near  0 
bearing  is  visible  only  near  the  start  and  end  of  the  run;  in  between  those  times  when  the 
scattering  patches  are  crossing  the  great  circle,  it  is  not  seen.  As  the  number  of  assumed  directions 
is  increased  from  5  to  7,  the  additional  directions  instead  of  indicating  the  sporadic-E  propagation 
are  appropriated  by  sidelobe  effects  similar  to  those  discussed  in  Section  6.2. 1 . 

The  sidelobes  of  the  10-element  mid-range  Vortex  array  are  relatively  large  at  the  8.0  MHz 
operating  frequency,  certainly  larger  than  those  of  the  Kestrel  array.  Thus  the  response  of  the  array 
when  steered  in  a  direction  having  a  sidelobe  directed  toward  a  scattering  patch  is  likely  to  be 
stronger  than  the  response  when  steered  toward  the  sporadic-E  direction.  In  addition,  the  narrower 
array  beamwidths  present  for  a  single  steered  direction  imply  that  a  scattering  patch  having  a 
particular  spread  in  direction  will  require  more  directions  to  be  allocated  to  it,  before  any 
remaining  direction  is  available  to  point  towards  a  weaker  sporadic-E  signal.  The  loss  of  visibility 
of  the  sporadic-E  mode  is  due  to  both  the  high  sidelobes  present,  and  the  narrower  beamwidths 
achieved  by  the  larger-aperture  10-element  Vortex  array. 

An  alternative  element  selection  was  made  from  the  Vortex  array,  to  correct  the  problem. 
The  innermost  eight  elements  were  selected.  This  resulted  in  a  spiral  radius  which  varied  from 
14.76  m  for  the  innermost  element  to  62.87m  for  the  outermost  element.  The  element  spacing 
covered  a  range  from  less  than  to  somewhat  greater  than  those  of  the  Kestrel  array.  The  sidelobes 
at  the  8  MHz  operating  frequency  could  be  expected  to  be  considerably  lower  than  the  previous 
element  selection. 

Figure  27  shows  the  MUSIC-estimated  bearing  and  elevation  time  histories  obtained  with 
the  8-element  Vortex  array  geometry,  for  5  assumed  signal  directions.  As  can  be  seen,  the 
sporadic-E  direction  is  very  well  defined  and  visible  throughout  the  run.  The  bearing  estimates  are 
tightly  clustered  about  the  true  bearing.  A  comparison  with  Figures  26  and  25  indicates  this  array 
configuration  to  be  the  best  for  the  conditions  simulated. 

The  results  of  this  section  show  that  sporadic-E  visibility  in  the  presence  of  stronger 
scattered  signals  is  highly  dependent  on  the  array  geometry.  Large  unfilled  arrays  with  high 
sidelobes  perform  poorly  in  this  regard.  Smaller  non-uniform  arrays  perform  very  well. 
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Figure  TJ.  Time  histories  of  estimated  bearing  and  elevation,  for  the  two-patch  geometry 
of  Figure  24  and  a  critical  frequency  of  7.0  MHz,  as  determined  by  the  MUSIC  algorithm 
using  the  Sinner  Vortex  array  elements,  for  5  assumed  signal  directions. 


6.2.4  Dependence  on  Relative  Signal  Strength 

One  factor  that  will  influence  how  often  sporadic-E  propagation  can  be  identified  in  the 
presence  of  larger  scattered  signals,  is  how  much  reduced  in  strength  the  sporadic-E  signal  can  be 
below  the  scattered  signals,  and  still  be  seen  by  the  MUSIC  algorithm.  In  order  to  determine  this, 
simulations  were  performed  for  the  single-scattering-patch  geometry  of  Figure  22,  using  the 
8-element  inner  Vortex  array  which  was  previously  identified  as  being  the  best  of  the  tested  array 
configurations.The  ionospheric  and  other  parameters  were  the  same  as  those  used  previously  in 
this  chapter,  with  the  exception  of  the  sporadic-E  reflection  coefficient  which  was  set  to  various 
levels  in  order  to  adjust  the  sporadic-E  signal  power  at  the  receiving  array  relative  to  the  scattered 
signal  power. 

Figure  28  shows  the  results  obtained,  in  terms  of  the  bearing  time  histories  for  relative 
power  levels  P(Es/scatter)  of  -11,  -15,  -21,  and  -26  dB. 

In  Figure  28,  when  its  power  is  -11  dB  relative  to  the  scattered  signals,  the  sporadic-E  trace 
is  extremely  well  defined,  and  correct  in  its  estimated  bearing.  When  its  power  is  reduced  further, 
to  -15,  -21,  and  finally  -26  dB  relative  to  scatter,  the  sporadic-E  becomes  less  well-defined  and 
precise  in  its  estimated  bearing.  At  -26  dB,  it  is  questionable  whether  the  sporadic  E  return  could 
be  identified  as  such,  especially  on  shorter  records  where  the  scattered  signals  may  not  show 
much  motion.  However,  at  -21  dB,  it  remains  sufficiently  clear  in  its  clustering  to  be  identified. 

Similar  tests  using  the  two-patch  geometry  of  Figure  24  show  that  the  sporadic-E 
propagation  yields  a  well-defined  estimated  bearing  trace  down  to  signal  levels  of  -17  dB. 
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P(Es/scatter)  =  -11  dB 


-15  dB 


Figure  28.  Estimated  bearing  time  histories  obtained  with  the  MUSIC  algorithm 
assuming  5  signal  directions,  for  the  single-patch  geometry  of  Figure  22  and 
various  sporc^ic-E  signal  powers  relative  to  the  scattered  signal. 


6.3  Summary  and  Discussion  of  Findings 

The  results  clearly  show  that  relatively  weak  sporadic-E  propagation  can  be  seen  in  the 
presence  of  stronger  scattering  propagation,  by  a  sampled-aperture  array  employing  the  MUSIC 
algorithm.  As  the  scattering  regions  provide  a  multiplicity  of  signal  paths  far  greater  than  the 
number  of  paths  that  can  be  assumed  by  MUSIC  with  a  small  number  of  array  elements,  the 
scattering  is  only  approximated  by  the  directions  found.  For  the  scattering  region  sizes  (200  km) 
and  the  array  apertures  considered  (50  m  at  8  MHz),  it  was  found  that  with  a  single  scattering 
region,  it  takes  at  least  two  assumed  directions  to  approximate  the  scattered  signals  sufficiently 
well  that  a  third  direction  could  provide  a  good  estimate  of  the  lower-power  sporadic-E  signal. 
Thus  with  a  single  scattering  patch,  MUSIC  must  assume  at  least  three  signal  directions,  and  with 
two  scattering  patches,  five  directions,  in  order  to  see  the  sporadic-E  return  weU  enough  to 
identify  it.  (Sporadic-E  can  be  identified  by  the  tightly  clustered  constant-bearing  estimates  and 
generally  low  estimated  elevations  characteristic  of  E-layer  propagation.) 

When  too  many  directions  are  assumed,  the  additional  directions  tend  to  be  appropriated  the 
sidelobes  of  the  array  response  to  the  stronger  signals.  This  is  especially  true  if  the  sidelobes  are 
close  in  amplitude  to  the  main  lobe.  The  resulting  false  traces  in  the  time  histories  create 
considerable  ambiguity,  making  sporadic-E  identification  difficult.  Also,  when  the  sidelobes  are 
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large  and  beamwidths  narrow,  as  for  the  larger-aperture  unfilled  10-element  Vortex  array  tested, 
the  estimated  directions  go  into  the  strong-scatter-signal  sidelobes  before  going  into  the  sporadic- 
E  direction.  The  sporadic-E  remains  undetected. 

These  results  present  implications  for  the  design  and  operation  of  DF  arrays.  While  laige- 
aperture  unfilled  DF  arrays  are  appropriate  for  separating  closely-spaced  signals  which  come 
from  separated  transmitters  on  one  or  several  propagation  paths,  their  narrow  beamwidths  and 
large  sidelobes  make  them  inappropriate  for  use  in  high-latitude  regions  where  extended 
scattering  regions  are  a  major  propagation  mechanism.  Smaller-aperture  arrays  for  which  signals 
scattered  from  an  extended  region  can  be  approximated  by  plane  waves  from  only  several 
directions  work  well  in  these  cases,  permitting  weaker  signals  such  as  sporadic-E  to  be  observed 
with  Uttle  ambiguity.  In  operating  such  arrays,  the  performance  depends  on  the  number  of 
assumed  signal  directions  used  by  the  DF  algorithm.  The  number  of  directions  required  to  account 
for  the  scattering  will  depend  on  the  number  and  extent  of  the  scattering  regions.  It  is  important 
for  the  operational  procedure  to  include  estimates  for  different  numbers  of  assumed  directions,  so 
that  the  optimum  for  any  particular  time  can  be  attained. 
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